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A.  SCIENTIFIC  OBJECTIVES:  The  objective  of  this  research  program  is  to  devise 
innovative  joint  time-frequency  (JTF)  processing  concepts  for  radar  image  enhancement 
and  physics-based  feature  extraction.  In  particular,  we  shall  investigate  how  JTF 
techniques  can  be  utilized  to  enhance  synthetic  aperture  radar  (SAR)  and  inverse 
synthetic  aperture  radar  (ISAR)  imageries  by  removing  artifacts  due  to  uncompensated 
target  motion,  complex  target  scattering  physics,  articulating  target  components  and 
clutter  and  propagation  effects.  Furthermore,  we  set  out  to  re-interpret  the  extracted 
artifacts  in  a  more  meaningful  feature  space  so  that  they  can  be  utilized  to  enhance  the 
performance  of  target  identification  algorithms.  This  research  is  leveraged  against  our 
previous  JTF  work  under  the  Joint  Services  Electronics  Program,  as  well  as  our  state-of- 
the-art  radar  signature  simulation  capabilities.  At  the  end  of  this  program,  we  plan  to 
develop  a  set  of  JTF-based  radar  image  processing  tools  in  the  standard  MATLAB 
environment,  so  that  this  work  can  be  easily  disseminated  to  the  radar  community. 

B.  SUMMARY  OF  RESULTS  AND  SIGNIFICANT  ACCOMPLISHMENTS: 

During  the  past  year,  we  have  proceeded  along  two  lines  of  research:  (i)  processing  real 
radar  data  to  test  our  JTF  algorithms  and  to  identify  needed  areas  of  research,  (ii) 
developing  new  algorithms  to  address  the  problem  areas  identified.  Along  the  first  line, 
we  have  continued  the  processing  of  measured  data  from  the  NATO  TIRA  radar.  By 


applying  the  JTF  motion  compensation  algorithm  we  had  developed  earlier,  we  have 
formed  high-resolution  ISAR  images  of  several  air  targets  and  have  compared  the  results 
to  truth  images  generated  using  instrumented  flight  data  and  simulated  images  computed 
using  the  radar  signature  prediction  code  Xpatch.  Through  these  comparisons,  we  have 
identified  several  problem  areas  for  blind  motion  compensation  due  to  the  presence  of 
either  non-steady  imaging  plane  or  moving  components  on  the  target.  These  issues 
cannot  be  handled  by  the  existing  motion  compensation  algorithms  and  are  being 
addressed  in  our  theoretical  algorithm  development.  We  have  also  initiated  an  effort  to 
process  the  NATO  MERIC  radar  data  set.  Blind  motion  compensation  is  more 
challenging  in  this  case,  as  the  flight  path  of  the  aircraft  is  very  close  to  the  radar,  which 
results  in  large  variations  in  the  aircraft’s  azimuth  and  elevation  angles  with  respect  to  the 
radar. 

Along  the  second  line  of  research,  we  have  focused  on  the  two  problem  areas 
identified  from  processing  of  the  TIRA  data.  First,  we  have  developed  and  applied  an 
adaptive  JTF  algorithm  to  remove  jet  engine  modulation  (JEM)  lines  caused  by  the 
scattering  from  the  rotating  engine  blades  of  the  aircraft.  The  algorithm  distinguishes  the 
slow  body  motion  from  the  fast  engine  blade  motion  using  the  phase  features  in  the  JTF 
space.  It  is  shown  that  JEM  lines  can  be  removed  from  the  image  to  unveil  the 
geometrical  features  of  the  aircraft.  Second,  we  have  initiated  research  to  address  the 
problem  of  blind  motion  compensation  in  the  presence  of  non-steady  imaging  plane. 
This  problem  can  be  considered  in  the  context  of  three-dimensional  motion  compensation 
and  is  much  more  challenging  than  conventional  two-dimensional  motion  compensation. 
Toward  this  end,  we  have  developed  an  algorithm  to  detect  the  presence  of  non-steady 
imaging  plane  during  the  imaging  interval  from  raw  radar  data.  This  detection  algorithm 
is  based  on  a  full  three-dimensional  motion  model  and  compares  the  phase  variation  of 
the  point  scatterers  in  different  range  cells  to  detect  the  presence  of  three-dimensional 
motion.  Furthermore,  we  have  begun  to  investigate  ways  to  form  a  focused  image  when 
the  target  undergoes  such  three-dimensional  motion.  We  believe  the  three-dimensional 
motion  compensation  theories  developed  in  this  topic  will  be  quite  important,  especially 
when  imaging  targets  that  exhibit  more  chaotic  motions,  such  as  ships  on  the  ocean 
surface.  The  detailed  descriptions  of  our  progress  are  described  below. 


Processing  of  NATO  TIRA  Radar  Data.  Our  research  thrust  during  the  first  year  of 
this  project  was  to  generate  ISAR  imagery  from  the  NATO  TIRA  radar  data  using  our 
previously  developed  adaptive  JTF  (AJTF)  algorithm  [1],  We  have  continued  this 
activity  into  the  second  year  by  comparing  the  results  against  truth  images  generated 
using  instrumented  flight  data  and  simulated  images  computed  using  the  radar  signature 
prediction  code  Xpatch.  While  the  comparison  among  the  images  is  quite  good,  we  have 
identified  several  needed  areas  of  research  in  ISAR-based  target  recognition.  Fig.  1 
shows  both  the  correlation  coefficient  between  the  JTF  and  truth  images  and  that  between 
the  JTF  and  Xpatch  images  versus  azimuth  look  angle.  From  the  two  curves,  we  see  that 
the  JTF  images  agree  very  well  with  the  truth  images,  indicating  that  blind  motion 
compensation  is  a  very  feasible  method  for  processing  real-world  radar  data.  The 
correlation  coefficient  between  the  JTF  and  Xpatch  images  is  slightly  lower  than  that 
between  the  measured  images.  In  particular,  two  problem  regions  can  be  seen  from  this 
plot.  First,  in  the  region  near  nose-on  (180  degrees  in  azimuth),  the  correlation  coefficient 
is  significantly  lower.  The  reason  is  due  to  JEM  lines  in  the  measured  data.  This  problem 
is  further  discussed  in  the  next  section  and  an  algorithm  to  remove  JEM  lines  is  proposed. 
Second,  at  some  angles  around  the  broadside  region  (90  degrees  in  azimuth),  the 
correlation  coefficient  is  also  low.  In  this  case,  both  the  associated  JTF  image  and  the 
truth  image  are  both  found  to  be  of  low  quality.  After  further  investigation,  it  is  found 
that  the  image  blurring  is  due  to  the  non-steadiness  of  the  imaging  plane.  This  problem  is 
also  being  further  addressed  in  our  research. 


Fig.  1.  Correlation  coefficients 
between  JTF  and  truth  images, 
and  between  JTF  and  Xpatch 
images. 


Removing  Image  Artifacts  due  to  Jet  Engine  Modulation.  Jet  engine  modulation 
(JEM)  is  a  well-known  phenomenon  caused  by  the  high-speed  rotation  of  the  aircraft 
engine.  For  imaging  radar,  the  typical  PRF  is  much  slower  than  the  engine  rotation 
frequency.  Thus  the  resulting  ISAR  image  in  the  frontal  region  of  an  aircraft  contains  an 
aliased  component  along  the  cross  range  dimension,  as  shown  by  the  TIRA  image  in  Fig. 
2(a).  Such  effect  is  difficult  to  predict  accurately  using  simulation  [2,  3].  Furthermore, 
JEM  lines  are  noise-like  and  can  corrupt  the  geometrical  features  of  the  target  in  the 
ISAR  image.  For  ISAR-based  target  recognition,  it  would  be  useful  to  devise  an 
algorithm  to  separate  the  JEM  lines  from  the  target  image  before  the  subsequent 
classification  process.  In  the  first  year  of  this  program,  we  successfully  demonstrated  the 
separation  of  the  rotating  blade  contribution  from  the  body  image  in  helicopter  data. 
However,  JEM  possesses  new  challenges  due  to  the  high  rotational  rate  of  the  jet  engine 
and  additional  electromagnetic  propagation  effect  through  the  inlet  duct.  We  have 
carried  out  JEM  removal  on  TIRA  data  using  the  AJTF  algorithm.  The  model  we  adopt 
assumes  that  the  aircraft  consists  of  a  slowly  rotating  body  with  a  constant  rotational 
velocity  Q*  and  a  fast  moving  engine  component  with  a  different  rotational  velocity  Qp. 
The  received  radar  return  as  a  function  of  dwell  time  can  thus  be  written  as: 


4/rf 
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where  N  is  the  total  number  of  point  scatterers  within  one  range  cell,  of  which  Nb  are  the 
body  scatterers.  Usually  Qp  is  much  greater  than  Qb-  While  the  first  term  can  be 
meaningfully  mapped  into  the  image  plane  of  the  target  via  the  Fourier  transform,  the 
second  term  results  in  serious  Doppler  smearing  across  the  cross  range  domain. 

We  can  also  utilize  the  AJTF  technique  to  separate  the  fast  moving  part  from  the 
relatively  slow  moving  body.  For  the  component  due  to  target  body  scattering,  the 
Doppler  frequency  is 

fb  _  W_ [ y cos(Q btD )  +  *  sin(fVD )]  =  —  nb (y  +  xQbtD ) 
c  c 

while  the  Doppler  frequency  due  to  the  fast  rotating  part  is 


(2) 
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Fig.  2.  TIRA  image  from  a  frontal  view,  (a)  Before  JEM  processing,  (b)  After 
JEM  removal  using  the  AJTF  algorithm. 

fS  =^-np[ycos(QptD)  +  xsm(nptD))  (3) 

We  see  that  (2)  is  a  linear  function  of  dwell  time  while  (3)  is  a  sinusoidal  function.  If  we 
parameterize  the  signal  by  basis  functions  that  have  linear  Doppler  frequency  behavior  as 
a  function  of  dwell  time,  the  two  signals  can  be  approximately  separated  by  their 
displacement  and  slope  parameters.  We  utilize  the  AJTF  processing  technique  to  carry 
out  the  parameterization.  The  signal  component  due  to  the  target  scattering  is 
reconstructed  by  using  all  the  bases  with  small  displacement  and  small  slope  parameters. 
Fig.  2(b)  shows  the  image  after  the  JEM  removal  processing.  Note  that  the  body  features 
are  unveiled  after  the  JEM  removal.  Fig.  3  shows  the  correlation  between  the  synthetic 
images  and  the  measured  images  from  TIRA  after  JEM  removal.  We  observe  that  the 
correlation  coefficients  in  the  JEM  region  are  increased  after  we  remove  the  JEM 
interference  from  the  body.  Further  research  is  needed  in  this  area  to  more  definitively 
assess  the  effect  of  JEM  removal.  This  will  be  carried  out  as  more  measured  data 
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Fig.  3.  Correlation  between  JTF  and  Xpatch  images 
after  JEM  removal. 


become  available.  Furthermore,  in  the  processing  of  the  MERIC  radar  data,  we  have 
observed  very  clear  indication  of  movement  from  a  mechanically  gimbaled  antenna  in  the 
nose  region  of  the  aircraft.  The  algorithm  will  be  adapted  to  study  this  issue. 

Three-Dimensional  Motion  Detection  Using  JTF  Algorithm.  One  basic  assumption  of 
existing  motion  compensation  algorithms  is  that  the  target  only  undergoes  motion  in  a 
two-dimensional  plane  during  the  dwell  interval  needed  to  form  an  image.  From  several 
independent  examinations  of  measured  ISAR  data  sets  recently,  it  was  reported  that  the 
presence  of  three-dimensional  motion  is  quite  detrimental  to  the  formation  of  a  well- 
focused  image  [5-7].  This  is  also  consistent  with  our  own  findings  from  the  TIRA  data. 
For  some  image  frames,  we  found  that  the  motion  compensated  image  and  the  associated 
truth  images  are  both  quite  poor.  An  example  image  is  shown  in  Fig.  4(a).  When  we 
examine  the  motion  data  from  the  aircraft  instrument,  we  find  that  the  aircraft  motion  is 
not  confined  to  a  two-dimensional  plane  during  the  imaging  interval.  Fig.  4(b)  shows  the 
plot  of  the  elevation  versus  azimuth  look  angles  of  the  aircraft  from  the  radar.  This  curve 
should  be  linear  if  the  motion  is  strictly  two-dimensional,  and  in  this  case,  there  is  clearly 
three-dimensional  motion.  In  general,  we  will  not  be  able  to  form  a  focused  image  using 
an  existing  motion  compensation  algorithm  based  on  the  two-dimensional  motion 
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Fig.  4.  Effect  of  non-steady  imaging  plane  on  image  quality,  (a)  A  poorly 
focused  image,  (b)  The  corresponding  target  motion  from  the 
instrumentation  data. 


assumption,  since  the  assumed  model  is  mismatched  to  the  actual  target  motion.  Our  goal 
in  this  research  is  to  develop  a  general  motion  compensation  algorithm  to  handle  targets 
with  arbitrary  three-dimensional  motion.  However,  this  problem  is  quite  challenging. 
Instead,  we  first  address  a  less  ambitious  problem  of  detecting  the  presence  of  three- 
dimensional  motion  from  raw  radar  data.  It  is  hoped  that  the  solution  here  will  be  a 
stepping  stone  to  the  ultimate  three-dimensional  motion  compensation  problem. 

Allowing  for  arbitrary  three-dimensional  motion  in  space,  we  consider  the 
following  model  as  a  generalization  of  the  two-dimensional  motion  model: 

N  4jjf 

£(t0J=!4  exp[-j-^H xk  +  yk6( tD  )  +  zk<p( tD  ))]  (4) 

k=i  c 

where  6  is  the  azimuth  angle  of  the  target  with  respect  to  the  radar,  and  (j)  is  the  elevation 
angle.  In  (4),  it  is  assumed  that  the  translation  motion  has  been  removed  and  that  the 
standard  small-angle,  small  bandwidth  approximations  apply.  This  model  reduces  to  the 
standard  two-dimensional  motion  model  when  0  and  (j)  are  linearly  related.  Our 
approach  to  the  three-dimensional  motion  detection  problem  is  to  utilize  the  AJTF 
algorithm  to  extract  the  phase  behavior  of  the  radar  data  at  multiple  range  cells.  It  can  be 
shown  that  when  the  target  undergoes  only  two-dimensional  motion  during  the  dwell 
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Fig.  5.  Blind  detection  of  three-dimensional  motion  from  TIRA  data,  (a) 
Degree  of  three-dimensional  motion  over  20  image  frames  detected 
using  the  proposed  algorithm,  (b)  Degree  of  three-dimensional 
motion  measured  from  the  instrument  data. 


duration,  the  relationship  between  the  phase  extracted  from  one  range  cell  and  that  from 
another  range  cell  should  be  linear.  For  three-dimensional  motion,  the  relationship  is  in 
general  nonlinear.  Therefore,  by  examining  the  linearity  of  the  phase  relationships  from 
different  range  cells,  we  can  distinguish  two-dimensional  motion  from  three-dimensional 
motion.  Fig.  5  shows  our  preliminary  results  from  the  TIRA  data.  Fig.  5(a)  shows  the 
degree  of  three-dimensional  motion  in  the  data  for  20  different  image  frames,  detected  by 
applying  our  algorithm  to  the  raw  TIRA  radar  data.  As  a  reference  for  comparison.  Fig. 
5(b)  shows  the  degree  of  three-dimensional  motion  for  the  same  20  frames  measured 


using  the  instrumentation  data.  It  can  be  seen  that  our  algorithm  correctly  detects  where 
significant  three-dimensional  motions  exist.  We  are  currently  fine  tuning  the  algorithm 
to  achieve  faster  and  more  robust  detection.  We  believe  this  detection  algorithm  could  be 
quite  useful  for  determining  the  “good”  imaging  intervals  from  which  focused  images  can 
be  more  readily  generated.  For  targets  that  exhibit  very  chaotic  motions,  such  as  ships  on 
the  ocean,  finding  such  intervals  of  opportunity  may  be  very  critical  for  target 
recognition.  Furthermore,  we  will  try  to  devise  algorithms  for  forming  focused  images 
even  in  the  presence  of  significant  three-dimensional  motions. 

Processing  of  NATO  MERIC  Radar  Data.  We  have  begun  to  process  the  MERIC 
radar  data  made  available  to  us  through  the  Naval  Research  Laboratory.  Blind  motion 
compensation  is  more  challenging  in  this  case,  as  the  flight  path  of  the  aircraft  is  very 
close  to  the  radar  and  the  instrument  data  show  more  variations  in  azimuth  and  elevation 
angles  during  the  flight  of  the  aircraft.  The  task  is  further  complicated  by  the  lack  of 
sufficient  instrument  data  (some  of  the  ARDS  pod  data  have  been  heavily  down-sampled 
and  written  into  Microsoft  Excel  form  in  the  MERIC  data  set).  In  addition,  some  of  the 
radar  data  contain  very  strong  colored  noise  in  the  background  and  considerable  effort 
has  been  spent  on  denoising  the  data.  Despite  these  difficulties,  the  final  JTF  image 
quality  we  have  generated  is  fairly  good.  Fig.  6  shows  the  angular  motion  parameters  of 
an  aircraft  for  which  the  radar  data  are  available.  The  aircraft  undergoes  an  azimuth 
rotation  of  about  25  degrees  during  the  17.5-sec  collection  period.  While  there  are 
35,000  pulses  of  radar  range  profiles  during  this  duration,  only  18  samples  of  instrument 
data  are  available  in  the  same  duration.  Consequently,  no  ground  truthing  was  possible 
due  to  the  lack  of  sufficient  instrument  data.  To  validate  our  results,  we  have  generated 
the  Xpatch  prediction  using  a  CAD  model  purchased  from  viewpoint.com.  It  is  however 
not  the  exact  model  for  the  aircraft  used  in  the  measurement.  The  model  was  first 
converted  into  an  Xpatch-compatible  format.  It  was  then  edited  to  remove  the  landing 
gear  and  bomb  loads  under  the  wings,  as  well  as  to  seal  up  the  open  seams  in  the  model. 
Some  preliminary  comparison  results  are  presented  in  Fig.  7.  Figs.  7(a)  and  7(b)  show 
respectively  the  JTF  motion  compensated  image  and  the  Xpatch  predicted  image  near 


point  P  of  Fig.  6  in  the  flight  path.  The  agreement  between  the  MERIC  image  and  the 
Xpatch  simulated  image  is  fairly  good,  considering  the  uncertainties  in  the  CAD  model. 


Figs.  7(c)  and  7(d)  show  respectively  the  JTF  motion  compensated  image  and  the  Xpatch 
predicted  image  near  point  Q  of  Fig.  6  in  the  flight  path.  It  is  clear  that  the  agreement  is 
quite  poor  in  this  case.  The  blind  image  formation  is  quite  challenging  during  this 
particular  portion  of  the  flight.  Furthermore,  we  are  not  certain  of  the  exact  imaging 
plane  for  the  simulation  due  to  the  sparseness  of  instrument  data.  We  plan  to  continue 
work  on  this  data  set  and  hope  to  eventually  use  the  data  set  to  further  our  investigation 
on  three-dimensional  motion  detection  and  compensation.  In  addition,  an  interesting 
phenomenon  was  observed  from  the  motion  compensated  image  of  a  second  aircraft.  The 
formed  images  show  Doppler  smearing  near  the  nose  region.  After  some  research  on  this 
aircraft,  we  have  determined  that  it  is  very  likely  due  to  movement  from  a  mechanically 
scanned  antenna  in  the  nose  cone  of  the  aircraft.  We  plan  to  examine  this  topic  in  more 
detail. 

C.  FOLLOW-UP  STATEMENT: 

During  the  coming  year,  our  research  effort  will  be  devoted  to  developing  new 
algorithms  to  attack  the  research  problems  we  have  identified  through  our  work  in 
measurement  data  processing.  The  research  areas  include:  (i)  three-dimensional  motion 
detection  and  three-dimensional  motion  compensation,  and  (ii)  image  formation  and 
feature  extraction  of  targets  with  multiple  moving  parts.  We  believe  these  two  problems 
are  of  fundamental  importance  in  ISAR  imaging,  as  the  current  state-of-the-art  ISAR 
algorithms  are  limited  by  the  two-dimensional  motion  and  rigid  body  assumptions.  For 
real-world  data  that  involve  more  complex  motions,  such  as  ships  with  arbitrary  motion 
and  moving  ocean  surface,  a  more  generalized  motion  compensation  framework  will  be 
highly  desirable.  In  addition  to  algorithm  development,  we  will  also  continue  our  effort 
to  process  new  radar  data  as  they  become  available  to  us.  Our  work  in  the  MERIC  data 
will  continue,  and  we  will  more  closely  examine  the  ship  data  from  the  Navy  small  craft 
ATR  program.  The  measurement  data  processing  enables  us  to  identify  relevant 
problems  that  can  feed  back  into  our  theoretical  investigation.  Furthermore,  the 
processed  data  sets  provide  us  with  a  well-understood  database,  from  which  we  can  test 
our  new  algorithms. 
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Adaptive  ISAR  Image  Construction  from  Nonuniformly 
Undersampled  Data 

Yuanxun  Wang  and  Hao  Ling 


Abstract — An  adaptive  approach  is  proposed  to  construct  ISAR  images 
from  nonuniformly  undersampled  data  in  the  angular  domain.  The  algo¬ 
rithm  uses  an  adaptive  scattering  feature  extraction  engine  in  place  of  the 
Fourier  transform  in  the  image  construction  procedure.  The  algorithm  en¬ 
tails  searching  and  extracting  out  individual  target  scattering  features  one 
at  a  time  in  an  iterative  fashion.  The  interference  between  different  target 
scattering  features  is  thus  avoided  and  a  dean  ISAR  image  without  the 
aliasing  effect  can  be  obtained.  The  algorithm  is  verified  by  constructing  the 
ISAR  image  from  the  chamber  measurement  data  of  the  model  VFY-218 
airplane. 

Index  Terms — Image  reconstruction,  synthetic  aperture  radar. 


I.  Introduction 

Constructing  an  inverse  synthetic  aperture  radar  (ISAR)  image  of  a 
target  requires  data  collection  in  both  the  frequency  and  angular  dimen¬ 
sions.  If  the  data  are  uniformly  sampled  and  the  sampling  rate  is  dense 
enough,  an  ISAR  image  can  be  obtained  by  using  a  two-dimensional 
(2-D)  fast  Fourier  transform  (FFT)  algorithm  [1].  In  this  paper,  we  ad¬ 
dress  the  case  when  the  angular  data  are  nonuniformly  undersampled. 
This  scenario  may  arise  in  real-world  data  collection  when  the  target 
is  fast  maneuvering  with  respect  to  the  radar  pulse  repetition  interval 
so  that  the  angular  look  on  the  target  by  the  radar  is  not  dense  enough 
to  satisfy  the  Nyquist  sampling  rate.  We  propose  an  algorithm  to  over¬ 
come  the  aliasing  effect  in  the  cross-range  dimension  and  construct 
ISAR  images  from  seriously  undersampled  data.  The  algorithm  uses 
an  adaptive  scattering  feature  extraction  engine  in  place  of  the  Fourier 
transform  in  the  image  construction  process.  The  original  concept  of 
adaptive  feature  extraction  was  proposed  in  [2]  and  [3].  It  has  been  ap¬ 
plied  to  ISAR  image  processing  in  the  joint  time-frequency  space  for 
resonant  scattering  mechanism  extraction  [4],  target  motion  compensa¬ 
tion  [5],  and  Doppler  interference  removal  [6].  In  contrast  to  the  Fourier 
transform,  where  the  signal  is  projected  to  all  the  image-domain  bases 
simultaneously,  the  adaptive  algorithm  searches  and  extracts  the  indi¬ 
vidual  target  scattering  features  one  at  a  time  in  an  iterative  fashion. 
When  applied  to  the  present  problem,  the  aliasing  error  caused  by  the 
interference  between  different  target  scattering  features  can  be  avoided. 
Therefore,  after  all  the  main  features  are  extracted,  they  can  be  dis¬ 
played  in  the  ISAR  image  plane  without  the  aliasing  effect.  We  verify 
this  algorithm  by  constructing  the  ISAR  image  using  the  chamber  mea¬ 
surement  data  of  the  model  VFY-218  airplane  [7].  It  is  found  that  a  rea¬ 
sonable  ISAR  image  can  be  constructed  from  seriously  undersampled 
data. 


II.  Adaptive  Feature  Extraction  (AFE)  Algorithm 

In  standard  ISAR  image  construction,  the  target  is  assumed  to  be  a 
collection  of  point  scattering  centers.  Under  the  small-angle  approx- 
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imation,  the  scattered  field  from  the  target  observed  as  a  function  of 
frequency  and  angle  can  be  written  as 

N 

E(f,8 )  =  £  0(xiyyi)e-2jkxi 
1  =  1 

=  £  0(xi,yi)e-2ikl^k^  (1) 

1=1 

where  0(x;,  y.)  is  the  amplitude  of  the  ith  scattering  center,  k  is  the 
free-space  wave  number,  and  kc  corresponds  to  the  wave  number  at 
the  center  frequency,  x,  and  y»  denote  the  down  range  and  cross-range 
dimensions,  respectively.  We  assume  that  the  sampling  in  frequency 
is  uniformly  spaced  and  dense  enough  to  satisfy  the  Nyquist  criterion 
since  it  is  completely  controlled  by  the  radar.  Thus,  the  range  pro¬ 
file  versus  angle  data  can  be  generated  from  the  frequency-aspect  data 
by  applying  a  one-dimensional  (1-D)  Fourier  transform  along  the  fre¬ 
quency  dimension.  We  shall  denote  the  result  as  R(x,  0) 

■  N 

R(x,  8)  =  J2  0(x„y,)Sz(x  -  X,)e-2>  W.  (2) 

1  =  1 

In  the  above  expression,  S*(x  -  x.)  is  the  down-range  point  spread 
function  due  to  the  finite-length  frequency  domain  data.  Similarly,  the 
cross-range  information  can  also  be  obtained  from  angular  data  via  a 
1-D  Fourier  transform  of  R(x ,  0)  along  the  angular  dimension.  The 
resulting  image  J(x,  y)  is 

I(x,y)  =  j  R{x,0)eiik'*e  d8 

N  . 

=  £  0(xi,yi)Sz(x  -  Xi)  I  e2]k^~^  d8 

N 

=  ^  ^  0{^Xi^yi)Sx{x  —  xi)Sy(y  —  yt)  (3a) 

i=i 

where 

Sv(y-yi)  =  J  e2ik^v~Vi)  d8  (3b) 

is  the  cross-range  point  spread  function  due  to  the  finite-length  angular 
domain  data.  If  the  data  are  sampled  densely  enough  such  that  the  nu¬ 
merical  integration  can  be  carried  out  accurately,  the  point-spread  func¬ 
tion  Sy  should  be  a  well-localized  function  with  its  peak  at  y,*,  while 
rapidly  decaying  away  from  the  peak.  The  resulting  image  J(x,  y)  will 
be  a  clean  image  with  good  indication  of  the  amplitudes  and  posi¬ 
tions  of  the  target  point  scattering  features.  However,  when  the  data 
are  undersampled,  the  numerical  integration  in  (3b)  will  result  in  large 
aliasing  error  that  shows  up  as  high  sidelobes  in  Sy.  Consequently,  the 
constructed  image  will  contain  strong  interference  between  the  scat¬ 
tering  features.  This  effect  can  be  interpreted  as  the  loss  of  orthogo¬ 
nality  of  the  Fourier  bases  under  the  undersampled  condition. 

In  the  proposed  approach,  we  use  an  adaptive  feature  extraction  al¬ 
gorithm  in  place  of  the  Fourier  processing.  Instead  of  projecting  the 
signal  onto  all  the  exponential  bases  simultaneously,  we  search  out  the 
strongest  point  scattering  feature  in  the  cross-range  domain  and  remove 
it  from  the  original  signal.  Then  the  search  is  repeated  for  the  remainder 
signal  and  the  point-scattering  features  are  extracted  one  at  a  time  until 
the  energy  of  the  residue  signal  is  smaller  than  a  preset  threshold.  The 
search  procedure  is  carried  out  by  calculating  the  integral  in  (3b)  for  all 
points  in  cross  range  but  saving  only  the  maximum  value  and  position, 
i.e.. 


[BP,yP]  =  max  [/P(x,y)J 

Vp 


(4) 


Fig.  1.  The  ISAR  image  constructed  at  30°  azimuth  from  the  original  data 
using  FFT. 


Fig.  2.  The  ISAR  image  constructed  at  30°  azimuth  from  randomly 
undersampled  data  using  Fourier  transform. 


where  p  denotes  the  stage  of  the  iterative  procedure.  The  remainder 
signal  is  produced  by  subtracting  out  the  pth  feature 

RP+l  (*,  8)  =  Rp(x,  6)  -  Bpt~2^9.  (5) 

The  convergence  of  the  above  procedure  is  guaranteed  and  the  math¬ 
ematical  proof  is  given  in  [2].  The  advantage  of  such  an  iterative  pro¬ 
cedure  is  that  each  time  we  extract  out  the  strongest  feature,  we  also 
eliminate  its  interference  on  the  other  features.  It  should  be  noted  that 
nonuniform  sampling  is  a  prerequisite  to  ensure  that  there  is  no  ambi¬ 
guity  in  the  strongest  features  since  uniformly  undersampled  data  will 
result  in  multiple  positions  of  the  strongest  features.  For  simplicity,  the 
algorithm  is  repeated  for  each  range  cell  of  the  image.  A  2-D  algorithm 
in  frequency  and  aspect  can  also  be  developed,  if  the  search  is  carried 
out  for  both  the  range  and  cross-range  parameters.  After  all  the  features 
are  extracted  out,  we  can  construct  an  ISAR  image  using  the  amplitudes 
and  positions  of  the  point  scatterers. 


III.  Results  and  Discussion 

To  examine  the  applicability  of  the  algorithm  on  real  target  scattering 
data,  we  reconstruct  the  radar  image  of  a  model  (1:30  scale)  VFY-218 
airplane  using  undersampled  chamber  measurement  data  [7].  The  mea¬ 
surement  data  consist  of  an  aspect  window  from  10°  to  50°  and  a  fre¬ 
quency  range  from  8  to  16  GHz.  To  construct  an  ISAR  image,  we  first 
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[7]  “Radar  cross  section  measurement  data  of  the  VFY  218  configuration,” 
Naval  Air  Warfare  Center,  China  Lake,  CA,  Tech.  Rep.  NAWCWPNS 
TM-7621,  Jan.  1994. 


Fig.  3.  The  ISAR  image  constructed  at  30°  azimuth  from  randomly 
undersampled  data  using  the  AFE  algorithm. 


polar  reformat  the  frequency-aspect  data  to  the  (Kx,  Ky)  space.  The 
reformatted  data  consist  of  401  samples  in  Kx  and  438  samples  in  Ky. 
The  ISAR  image  is  first  generated  by  fast  Fourier  transform  (FFT)  and 
shown  with  the  airplane  overlay  in  Fig.  1.  The  point-scattering  features 
can  clearly  be  seen.  Next,  we  test  our  algorithm  by  generating  an  under¬ 
sampled  data  set  in  Kv.  This  is  approximately  the  same  as  undersam¬ 
pling  in  angle.  (Note  that  for  full  size  targets,  this  approximation  gets 
better.)  The  Nyquist  sampling  rate  requires  about  36  sampling  points 
in  Ky  and  we  randomly  select  only  24  out  of  the  438  points.  The  max¬ 
imum  sampling  interval  used  is  about  four  times  the  size  of  the  Nyquist 
sampling  interval.  Therefore,  serious  aliasing  will  occur  if  the  Fourier 
transform  algorithm  is  used,  as  shown  by  the  ISAR  image  displayed  in 
Fig.  2.  All  the  features  are  overlapped  with  sidelobe  noise  such  that  the 
point  scattering  features  on  the  airplane  can  no  longer  be  distinguished. 
Next,  we  apply  the  adaptive  feature  extraction  (AFE)  algorithm  to  each 
range  cell  of  the  image.  The  image  is  reconstructed  and  shown  in  Fig.  3. 
Comparing  Fig.  3  with  Fig.  1,  we  can  see  the  main  features  of  the  air¬ 
plane  in  Fig.  3  are  all  well  reconstructed  in  Fig.  1 .  We  do  observe  some 
noisy  spots  outside  the  target  at  the  lower  dynamic  ranges  in  Fig.  3. 
This  low-level  noise  occurs  at  about  25  dB  down  from  the  key  features 
and  presents  a  dynamic  range  limitation  of  the  present  AFE  algorithm. 
The  algorithm  has  also  been  tested  on  measured  data  from  in-flight  tar¬ 
gets  with  good  success. 
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Abstract-An  adaptive  wavelet  packet  transform  (AWPT)  algorithm 
is  proposed  to  process  synthetic  aperture  radar  (SAR)  imagery  and 
remove  background  clutters  from  target  images.  Since  target  features 
are  more  efficiently  represented  using  the  wavelet  packet  bases,  higher 
signal-to-clutter  ratios  (SCR)  can  be  achieved  in  the  wavelet  transform 
domain.  Consequently,  clutters  can  be  more  effectively  separated  from 
the  desired  target  features  in  the  transform  domain  than  in  the  orig¬ 
inal  SAR  domain.  The  processed  results  based  on  the  MSTAR  data 
set  demonstrate  the  effectiveness  of  this  algorithm  for  SAR  clutter  re¬ 
duction. 
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1.  INTRODUCTION 

Synthetic  aperture  radar  (SAR)  images  of  ground  targets  after  the 
focusing  process  generally  consist  of  target  images  and  clutters  from 
background  scattering  [1].  The  clutter  signal  is  undesirable  since  it  in¬ 
terferes  with  the  actual  target  features  in  automatic  target  recognition 
(ATR)  applications,  and  has  to  be  removed  before  ATR  processing. 
The  traditional  way  to  suppress  clutter  is  to  choose  a  threshold  using 
CFAR  that  is  based  on  the  standard  deviation  and  mean  value  of  the 
clutter,  and  the  threshold  is  applied  to  the  SAR  image  chip  to  remove 
the  clutter  directly  [2,  3].  However,  this  approach  assumes  that  the 
target  signal-to-clutter  ratio  (SCR)  is  large  enough.  If  this  assumption 
is  not  true,  this  approach  results  in  either  target  feature  loss  or  large 
clutter  residues.  In  this  work,  we  investigate  the  use  of  orthogonal  basis 
transformation  to  increase  the  SCR  for  more  effective  clutter  removal. 
More  specifically,  our  approach  is  to  find  the  best  wavelet  packet  basis 
to  represent  a  SAR  image  using  the  adaptive  wavelet  packet  transform 
(AWPT)  [4-9]. 

Wavelet  packet  basis  is  a  generalization  of  the  conventional  wavelet 
basis.  It  retains  the  multi-resolution  property  of  the  conventional 
wavelet  basis  while  abandoning  the  rigid  constant-Q  structure  in  the 
decomposition.  Wavelet  packets  have  been  applied  to  image  compres¬ 
sion  [7]  and  moment  matrix  sparsification  [9]  with  good  success.  Theo¬ 
retically  they  include  regular  pulse  basis,  FFT  basis,  short-time  Fourier 
transform  (STFT)  basis,  as  well  as  conventional  wavelet  basis.  In  this 
application,  we  transform  the  SAR  image  into  the  wavelet  packet  basis 
to  maximize  the  signal-clutter-ratio  in  the  transform  domain.  Since 
a  typical  target  image  usually  consists  of  point  scatterers  and  more 
diffused  region  features,  the  multi-scaled  wavelet  basis  is  well  suited  to 
focus  the  target  image.  Clutter  in  the  images,  on  the  other  hand,  is  sta¬ 
tistically  uncorrelated  (or  weakly  correlated)  from  pixel  to  pixel,  and 
the  transformed  clutter  image  under  the  same  set  of  bases  remains  un¬ 
focused  [8, 10, 11].  Therefore,  we  expect  that  the  SCR  can  be  increased 
by  transforming  the  original  image  with  an  appropriately  chosen  set  of 
wavelet  packet  basis.  To  determine  the  best  basis  for  a  SAR  image,  we 
implement  the  AWPT  algorithm  to  search  for  the  best  wavelet  packet 
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basis  based  on  a  cost  function  that  describes  how  well  the  target  signal 
is  focused  in  the  transform  domain.  Because  clutter  in  a  SAR  image  is 
usually  not  strictly  white,  we  apply  a  frequency-dependent  threshold 
to  the  transformed  image.  We  then  inverse-transform  the  thresheld 
image  in  the  wavelet  packet  basis  domain  back  to  the  SAR  domain  to 
obtain  a  de-cluttered  target  image. 

This  paper  is  organized  as  follows.  In  Section  2,  we  introduce  the 
basic  concept  of  wavelet  packet  basis.  In  Section  3,  we  discuss  problem 
.formulation  and  the  algorithm  for  SAR  clutter  rejection  based  on  the 
adaptive  wavelet  packet  transform.  In  Section  4,  we  present  the  test 
results  of  the  algorithm  using  the  MSTAR  data  set  [12].  In  Section  5, 
we  draw  some  conclusions  from  this  work. 

2.  WAVELET  PACKET  BASIS  FOR  SIGNAL 
REPRESENTATION 

A  wavelet  packet  basis  function  can  be  expressed  in  the  space  domain 
as: 

<Pjtk(x)  =  2-j/Vn(2-Jx -k),  he  z,  jez,  n  e  z+  (l) 

where  k,  j,  and  n  denote  the  space  shift,  scale,  and  modulation  index, 
respectively.  The  function  if>n  can  be  generated  from  the  decomposi¬ 
tions  of  both  the  scaling  function  and  mother  wavelet  function  using 
the  “two-scale  equation”  [4,  5].  The  parameter  choice  of  a  wavelet 
packet  basis  is  not  unique  for  a  complete  and  orthogonal  decomposi¬ 
tion  of  a  signal  in  L2  space.  For  a  wavelet  packet  basis  function 
we  define  its  dyadic  interval  /J>n  C  3?  in  the  frequency  domain  as  [4, 
9]: 

/;,«=[ 2-'n,  2"J(n  +  l))  (2) 

For  the  complete  and  orthogonal  wavelet  basis  functions  {<£>”}  chosen 
to  represent  a  signal  in  L 2  space,  their  dyadic  intervals  should  be 
disjoint  and  cover  the  entire  signal  bandwidth,  i.e., 


Pl-O.n  —  [Q)\ 

3,n 

(3) 

U4n  =  (0,l) 

j,n 

(4) 

and 
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Therefore,  the  wavelet  packet  basis  function  <pf  can  be  interpreted  as 
a  modulated  wavelet  with  a  central  frequency  of  2 -J(n  4-  1/2)  and  a 
normalized  bandwidth  of  2-J. 

For  two-dimension  image  decomposition,  the  2-D  wavelet  packet 
basis  function  can  be  configured  from  two  1-D  wavelet  packet  basis 
functions  as  follows: 


Uiq(m,n)  =  (5) 

In  (5),  we  restrict  the  scales  of  the  two  1-D  basis  functions  to  be  the 
same.  This  guarantees  the  basis  functions  will  be  of  the  same  size  in 
the  horizontal  and  vertical  directions,  and  remove  the  inter-scale  cou¬ 
pling  terms  in  the  transformed  image.  For  a  complete  and  orthogonal 
decomposition,  the  two  1-D  basis  functions  in  (5)  must  also  meet  the 
requirements  in  (3)  and  (4). 

Assuming  that  the  original  image  in  the  spatial  (SAR)  domain  is 
{s(m,  n),  0  <  m,  n  <  N},  we  define  a  set  of  complete  and  orthonormal 
2-D  wavelet  packet  basis  functions: 

{Uiq(k,  l)  0  <  j  <  J,  0  <  p,  q  <  2’,  0  <  k,  l  <  AT2-*}  (6) 

where  j  denotes  the  scale,  J  =  log2(iV),  p  and  q  are  the  frequency 
modulation  indices,  and  k,  l  are  the  position  indices.  The  decompo¬ 
sition  coefficients  of  the  image  s(m,  n)  using  the  wavelet  packet  basis 
(6)  are: 


sL(k>  o  =  E  E  nWfc  -  2im> 1  - 2in)  w 

m  n 

With  the  image  extended  to  be  periodic  in  the  spatial  domain,  the 
total  number  of  coefficients  is  exactly  N2  through  the  transform  in 
(7).  For  convenience,  we  refer  to  the  coefficients  as  simply  a  matrix 
{5(m,  n),  1  <  m,  n  <  N }  without  explicit  indication  of  j,  p  and  q. 

Given  a  SAR  image  s(m,  n),  we  set  out  to  find  the  best  wavelet 
packet  basis  to  maximize  the  SCR  in  transform  domain.  In  the  other 
words,  we  try  to  make  the  transformed  coefficient  matrix  {S(m,n), 
1  <  m,  n  <  N}  as  sparse  as  possible. 
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3.  SAR  IMAGE  CLUTTER  REDUCTION  USING 
ADAPTIVE  WAVELET  PACKET  TRANSFORM 


3.1  Statistical  Signal  Model 

We  assume  that  a  focused  SAR  image  in  spatial  domain  is: 


{s(m,  n),  0  <m,n<  N} 


(8) 


Ignoring  additive  noise  and  possible  multiplicative  speckle  noise  in  the 
image,  we  consider  a  SAR  image  that  consists  of  only  the  target  image 
and  clutters  generated  from  the  background  scattering  [15].  Basically 
the  target  image  and  the  clutters  are  non-overlapping  in  a  SAR  image, 
therefore  there  is: 


f  t(m,  n),  (m,  n)  in  target  area 

\  c(m,  n),  (m,  n)  in  clutter  area 


0  <  m,  n  <  N 


(9) 


For  a  clutter  pixel  c(m,  n),  it  includes  reflections  from  minute  scatters 
inside  this  resolution  cell.  For  modem  coherent  imaging  radar  with  a 
quadrature  receiver,  it  can  be  represented  as: 


K 

c(m,  n)  =  ^2akexp(j2nfcTk) 
k= l 

=  Cr  +  jCi 


(10) 


where  ak  and  rk  are  the  magnitude  and  relative  two-way  propagation 
delay  of  the  reflection  from  the  fc-th  scatter  in  this  clutter  resolution 
cell.  If  we  assume  those  variables  are  statistically  independent  and  K 
is  large,  according  to  Central  limit  theorem  the  sum  in  (10)  tends  to 
be  a  complex  Gaussian  random  variable,  and  its  amplitude  is  approxi¬ 
mately  Rayleigh  distributed  [13,  16,  17].  Although,  for  high-resolution 
radar  with  a  low  gazing  angle,  the  amplitude  is  closer  to  be  Weibull 
distributed  for  ground  clutters,  Rayleigh  approximation  is  accurate 
enough,  especially  for  the  distribution  near  the  mean  [15]. 

Another  important  assumption  about  SAR  clutters  is  that  clutters 
in  different  pixels  are  statistically  independent  [13-16].  Therefore,  we 
have: 


E  [c(m,  n)c*(p,  ?)]  =  a%6(p  -  m)5(q  -  n) 


(11) 
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where  ac  is  the  clutter  standard  deviation,  and  S  is  the  Dirac  function. 

Using  ergodic  characteristics  of  independent  identical  random  vari¬ 
ables,  the  statistical  mean  value  in  (11)  can  be  approximated  with  a 
spatial  average: 

E  [c(m,  n)c*(p,  ?)]  -  0c*(fc  +  P-mJ  +  Q-n) 

k  l 

=  a^6(p-m)6(q-n)  (12) 

From  (12)  we  find  the  clutters  are  spatially  uncorrelated,  which  provide 
clutters  a  distinct  feature  from  the  target  signal.  But  we  must  note 
that  for  SAR  clutters  there  is  some  weak  correlation  among  adjacent 
pixels,  thus  (12)  is  not  quite  right  if  p  —  m  and  q  —  n  are  close  to 
zeros. 

Generally  the  target  area  is  in  the  central  part  of  a  SAR  image,  and 
the  clutters  are  located  in  surrounding  parts.  We  define  that  the  target 
signal  is  zero  in  clutter  area,  and  the  clutter  is  zero  in  target  area,  i.e., 

t(m,  n)  =  0,  (m,  n)  in  clutter  area  (13) 

c(m,  n)  =  0,  (m,  n)  in  target  area  (14) 

Because  zero  values  don’t  contribute  anything  to  the  summation  in 
(12),  the  definition  in  (14)  doesn’t  affect  the  autocorrelation  property 
of  clutters  in  (12)  as  long  as  there  axe  enough  actual  clutter  samples  in 
the  correlation.  Also  the  zeros  in  (13)  won’t  affect  the  autocorrelation 
property  of  target  image. 

Considering  (9),  (13),  and  (14),  we  have: 

s(m,  n)  =  t(m,  n)  +  c(m,  n)  0  <  m,  n  <  N  (15) 

With  the  additive  model  in  (15)  for  a  SAR  image,  we  need  to  find  the 
best  wavelet  packet  basis  to  maximize  Signal-to-Clutter  ratio  in  the 
transform  domain. 

3.2  Best  Wavelet  Packet  Basis  for  Maximization  of  Signal-to- 
Clutter  Ratio 

If  we  apply  a  wavelet  packet  basis  shown  in  (6)  to  a  SAR  image 
s(m,n)  in  (15),  the  transform  coefficients  are: 

sL(k>  o  =  EE  n)uUk  -  2'm’ 1  -  2'n) 

m  n 

=Ti,(k,l)  +  Cl,{k,l) 


(16) 
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where: 


!?,,(*,  o  =  E£  *(•»■  nH«(k  - 2>m ■ '  - 2J")  <17) 

m  n 

% ,(k,l)  =  ££>(».,  n)c£,(fc  -  2i">  <18) 

m  n 

are  the  coefficients  for  target  image  and  clutters,  respectively. 

Because  a  liner  combination  of  independent  Gaussian  random  vari¬ 
able,  the  clutter  coefficients  still  are  Gaussian  distributed.  If  we  re¬ 
write  (18)  into  a  simple  form,  and  note  that  wavelet  packet  basis  is 
orthonormal,  there  are: 

C(k,  l )  =  J2  c(m’ n)^  (k»  h  m) n)  (19) 

m  n 


and 


s[c(fc,o^(i,i)] 


=  £? 


c(m,  n)C7 (fc,  Z,  m,  n)  ^  ^  c*(p,  g)Z7(i,  j,  p,  q) 


V  Q 


=  ££££  E  [c(m,  n)c*(p,  g)]  t/(fc,  Z,  m,  n)U(i,j,p,  q) 

m  n  p  q 

=  ££££  ac6(p  -  m)5(q  —  n)U(k,l,m,n)U(i,j,p,q) 

m  n  p  q 

Z,  m,  n)U(i,  j ,  m,  n) 

m  n 

=  <rc<S(t  -  fc)«J(i  -  0 


(20) 


Therefore,  the  clutter  coefficients  are  uncorrelated,  thus  independent 
Gaussian  statistical  variables  with  the  same  mean  and  variance  as  those 
of  the  clutter  in  spatial  domain. 

We  have  demonstrated  that  the  statistical  characteristics  of  SAR 
image  clutters  don’t  change  after  wavelet  packet  basis  transformation. 
But  target  image  is  correlated  to  itself.  We  need  to  find  a  wavelet 
packet  basis  to  make  transformed  target  image  further  concentrated, 
and  increase  the  signal-to-clutter  ratio  in  transform  domain,  in  other 
words,  we  need  to  make  transformed  target  image  sparse.  A  cost  func¬ 
tions  is  needed  to  measure  the  sparsity  of  the  transformed  signal,  thus 
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the  best  wavelet  packet  basis  can  be  found  for  a  specific  image  to 
achieve  the  maximum  sparsity  in  the  transformed  image.  A  typical 
cost  function  for  sparsity  is  entropy  function.  But  entropy  function 
itself  is  not  an  additive  cost  function,  which  makes  it  inapplicable  in 
a  fast  global  best  basis  search  algorithm.  A  modified  additive  en¬ 
tropy  function  was  proposed  [18,  19],  but  it  requires  computation¬ 
intensive  logarithm  operation  to  evaluate  function  cost,  and  is  unsta¬ 
ble  in  the  environment  with  clutter  or  noise  interference.  Alternatively 
we  use  lp  energy  concentration  function  as  cost  function  in  this  ap¬ 
plication.  The  lp  energy  concentration  function  of  a  data  sequence 
(x(&),  k  =  1,2, .  ..,AT}  is  defined  as: 

K 

C  =  J2\x(k)\*  0  <p  <2  (21) 

k=  1 

For  a  transformed  SAR  image  with  a  coefficient  matrix  ,  the  cost 
function  is:  _ 

C  =  EE  0  <  p  <  2  (22) 

k  l 

For  convenience  and  efficiency,  we  usually  choose:  p  =  1  in  (22).  Sup¬ 
posing  we  transform  a  sequence  {xi,  X2,  •  •  • ,  zl)  into  another  sequence 
{xi,  X2,  ■  •  • ,  ~xL}  with  a  wavelet  packet  basis  using  the  cost  function  of 
(22),  the  transformed  sequence  is  the  sparsest  if  all  of  them  excerpt 
one  are  zeros,  accordingly  the  cost  is  the  minimum  in  this  case.  We 
illustrate  the  best  basis  selection  criterion  based  on  the  energy  concen¬ 
tration  cost  function  in  (21)  using  a  simple  case:  assuming  that  a  data 
sequence  is  {xi,X2}  with  energy  equal  to  1,  and  its  orthonormal  basis 
transform  coefficients  are  {xi,X2}  ,  we  have: 

N2  +  N2  =  l  (23) 


The  cost  function  is: 


Cost  =  |xi|  +  |x2|  (24) 

As  shown  in  Figure  1,  transformed  data  must  be  a  point  on  the  circle 
with  radius  of  1  because  of  constant  energy  constraint.  The  Lines 
u,  b,  and  c  are  the  energy  concentration  cost  functions  with  p  =  1. 
The  cost  function  lines  must  cross  the  one-fourth  of  the  circle  to  satisfy 
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Figure  1.  Relationship  between  the  sparsity  of  a  basis  transformed 
signal  and  the  energy  concentration  cost  function. 

the  energy  constraint.  Apparently  on  Line  c  the  blank  point  is  the 
transformed  data,  but  it  has  the  highest  cost,  and  worst  sparsity.  On 
Line  a,  two  white  points  axe  the  possible  transformed  data.  They  have 
the  lowest  cost,  and  best  sparsity  for  all  possible  basis  transformation. 
The  gray  points  on  Line  b  are  the  transformed  data  with  cost  and 
sparisty  between  the  best  and  the  worst.  Therefore,  based  on  the 
energy  concentration  cost  function,  we  can  effectively  find  the  best 
transform  basis  to  maximize  the  sparsity  of  the  transformed  signals. 

Due  to  the  non-overlapping  of  target  and  clutter  in  spatial  domain, 
mostly  the  coefficients  for  target  and  clutter  are  also  separated  in  the 
finer  scales  in  the  transform  domain.  Thus  the  total  cost  of  the  trans¬ 
formed  SAR  images  in  (16)  is  approximately  equal  to  the  sum  of  the 
costs  of  transformed  target  images  and  clutters,  i.e., 

Cost  ( [§] )  a*  Cost  (  [t]  )  +  Cost  ( [6] )  (25) 

where  sj  ,  [VJ ,  and  |cj  are  the  transform  coefficient  matrices 
for  the  whole  SAR  image,  target  image,  and  clutters,  respectively. 
The  cost  of  transformed  clutters  doesn’t  change  because  they  are  the 
random  variables  with  exactly  the  same  distribution  as  that  of  the 
original  clutters  in  spatial  domain.  The  minimization  of  the  cost  of 
the  whole  transformed  SAR  image  is  equivalent  to  the  minimization 
of  that  of  the  transformed  target  image,  and  thus  the  maximization  of 
target  signal-to-clutter  ratio  in  the  transform  domain. 
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3.3  Implementation  of  Adaptive  Wavelet  Packet  Transform 
for  SAR  Images 

Theoretically  we  can  apply  every  possible  wavelet  packet  basis  to  a 
SAR  image,  evaluate  the  cost  of  the  every  transformed  image,  and  find 
the  best  wavelet  packet  basis  with  the  least  cost.  But  the  brutal  force 
computation  is  impractical  to  implement  because  there  are  too  many 
eligible  wavelet  packet  bases  for  a  specific  application  and  the  direct 
convolution  in  (7)  is  computationally  too  costly.  Fortunately,  in  [18]  a 
fast  decomposition  algorithm  was  proposed  to  search  the  best  wavelet 
packet  for  a  data  sequence  based  on  additive  cost  function,  and  it  can 
be  readily  adapted  to  2-D  SAR  image  processing. 

Generally  the  original  SAR  image  is  the  discretized  samples  using 
pulse  basis,  and  it  can  be  approximated  as  the  coefficients  of  wavelet 
packet  basis  with  the  highest  spatial  resolution,  i.e.,  in  the  finest  scale. 
The  decomposition  coefficients  in  the  next  scale  are  related  to  the  ones 
in  current  scale  by  “2-scale  equations”  using  a  pair  of  quadrature  filters 
{h(n)}  and  {^(n)}. 

Assuming  that  the  initial  image  samples  are  represented  with  a  ma¬ 
trix  [S],  its  decomposition  coefficients  at  the  next  scale  can  be  obtained 
through  the  convolution  and  down-sampling  of  [S]  and  quadrature  fil¬ 
ter  impulse  responses  {/i(n)}  and  (g(n)},  i.e., 

Stf  Jr(m ,  n)  =  ^  s(k ,  l)g(2m  -  k)g(2n  - 1 ) 

k  l 

n )  =  EE  s(k,  l)h{2m  —  k)g(2n  —  l) 

k  1  (26) 
SHL(m ■> n)  ~  XI ~  k)h(2n  -  l ) 
fc  l 

EE  s(k,l)h(2m-k)h(2n-l) 

k  l 

where  {h(n)}  and  {g(n)}  denote  the  impulse  responses  of  low-  and 
high-pass  Quadrature  Filters  (QFs),  respectively.  The  decomposition 
in  (26)  is  widely  referred  as  quadtree  decomposition,  which  is  illus¬ 
trated  in  Figure  2.  The  quadtree  decompositions  can  be  applied  recur¬ 
sively  until  the  coarsest  scale  is  reached.  Apparently  the  decomposition 
at  each  scale  needs  0(N2)  operations,  and  the  number  of  total  avail¬ 
able  scales  is  about  log2(N).  Therefore,  the  total  computation  cost 
to  implement  the  full  decomposition  from  spatial  domain  to  frequency 


Clutter  reduction  for  synthetic  aperture  radar  imagery 


11 


S% 

C<2) 

*2 

C<2) 

“JOT 

Figure  2.  Quadtree  decomposition  structure. 


domain  is  about  0(N2\og2N)  operations  using  the  quadtree  decompo¬ 
sition  method.  Because  the  cost  function  is  additive,  the  best  wavelet 
packet  basis  and  corresponding  transform  coefficients  can  efficiently 
found  from  the  fully  quadtree  decomposed  result.  With  regard  to  a 
SAR  image  the  procedures  to  find  the  best  wavelet  packet  basis  and 
the  transform  coefficient  are  as  follows: 

1.  Fully  decompose  the  initial  image  into  frequency  domain  using 
quadtree  decomposition. 

2.  Use  the  wavelet  basis  at  the  last  scale  as  the  best  initial  basis, 
and  its  cost  as  the  best  initial  cost. 

3.  Trace  back  from  the  last  scale,  and  comparing  the  cost  at  every 
node  with  the  smallest  cost  in  the  all  branches  decomposed  from 
the  node.  If  the  cost  is  reduced,  update  the  best  basis  and  the 
corresponding  cost  using  this  node;  otherwise  continue  the  back¬ 
ward  search  for  the  best  basis. 

4.  The  best  basis  and  its  transformed  coefficients  axe  found  when 
the  search  comes  toward  the  first  scale,  i.e.,  the  initial  spatial 
sampling  data. 


4.  PROCESSING  RESULTS 

With  the  principles  described,  the  clutter  in  a  specific  SAR  image  can 
be  removed  in  the  wavelet  packet  basis  domain.  The  clutter  removal  is 
implemented  by  thresholding  processing  in  the  transform  domain  with 
a  higher  Signal-to-Clutter  ratio  compared  with  direct  thresholding  in 
spatial  domain.  But  even  the  thresholding  processing  in  transform 
domain  will  inevitably  cause  the  target  image  loss,  while  the  clutter 
is  suppressed.  In  the  processing,  we  choose  a  threshold  level  to  keep 
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target  image  loss  under  an  acceptable  level,  and  allow  some  clutter 
residues  to  exist  after  the  processing.  The  clutter  residues  can  be 
easily  removed  by  some  post-processing  such  as  additional  thresholding 
and  clustering  processing  in  spatial  domain.  The  processing  scheme  for 
SAR  clutter  reduction  using  AWPT  method  is  shown  in  Figure  3.  This 
clutter  reduction  scheme  is  applied  to  standard  MSTAR  SAR  images 
[12].  A  typical  focused  MSTAR  image  with  clutters  is  shown  in  Figure 
4.  The  image  size  is  128  x  128,  and  the  resolution  is  about  0.3  m  x 
.0.3  m.  The  target  in  the  image  is  BTR-70  transportation  vehicle,  and 
the  clutter  was  generated  from  background  vegetation  scattering.  From 
the  image,  we  can  find  that  the  reflections  from  the  rear  part  of  this 
vehicle  are  weak  and  indistinguishable  from  the  background  clutters 
due  to  the  signal  blockage  by  the  front  part  of  the  target.  It  would  be 
difficult  to  remove  all  clutters  and  keep  the  whole  target  signal  intact 
using  direct  thresholding  method. 

4.1  Wavelet  Filter  for  SAR  Image  Transform 

Wavelet  packet  transform  of  an  image  can  be  implemented  by  re¬ 
cursive  quad-tree  decomposition  of  the  original  spatial  sampling  data. 
Based  on  “two-scale  equations,”  the  coefficients  in  the  coarser  scales 
can  be  derived  from  the  coefficients  in  the  finer  scales.  Considering 
the  original  sampling  SAR  data  as  the  transform  coefficients  with  the 
finest  scale  in  spatial  domain,  we  can  obtain  all  wavelet  packet  trans¬ 
form  coefficients  by  repeatedly  filtering  and  down-sampling  of  the  orig¬ 
inal  spatial  image  using  a  wavelet  filter.  To  maximally  concentrate  the 
energy  of  the  transformed  image,  we  choose  Daubechies  wavelet  filter 
in  this  application.  It  is  the  most  commonly  used  orthonormal  wavelet 
filer;  and  the  corresponding  wavelet  basis  functions  possess  localiza¬ 
tion  property  in  both  time  and  frequency  domains.  The  order  of  the 
wavelet  filter  is  related  to  the  vanishing  moments  of  the  wavelet  ba¬ 
sis  functions.  The  higher  the  wavelet  filter  order,  the  more  vanishing 
moments  the  basis  functions,  thus  the  more  concentrated  the  trans¬ 
formed  coefficients.  However,  through  the  convolution  operation,  the 
finest  spatial  resolution  available  for  the  transformed  coefficients  is  the 
order  of  the  wavelet  filter.  For  a  SAR  image  with  scattering  points 
in  high  spatial  resolution,  the  wavelet  filter  of  low  order  is  desired  for 
efficient  representation  of  those  scattering  points  in  the  transform  do¬ 
main.  Therefore,  we  choose  Daubechies  filter  with  the  order  of  6  as  the 
wavelet  filters  in  this  application.  Using  the  AWPT  algorithm  based  on 
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Figure  3.  The  processing  scheme  of  automatic  SAR  clutter  reduction 
based  on  Adaptive  Wavelet  Packet  Transform  (AWPT). 


Figure  4.  SAR  image  for  BTR-70  transportation  vehicle  with  vegeta¬ 
tion  clutters. 

the  quad- tree  decomposition,  we  transformed  BTR-70  image  in  Figure 
4  into  AWPT  domain.  The  transformed  image  is  displayed  in  Figure 
5.  The  value  of  the  corresponding  energy  concentration  cost  function 
is  about  595  for  the  best  wavelet  packet  basis,  opposed  to  777  of  the 
cost  function  value  in  spatial  domain;  while  the  cost  function  value  is 
about  641  if  the  conventional  wavelet  transform  is  applied  to  the  im¬ 
age.  Apparently  we  achieve  the  sparsest  transformed  image  using  the 
AWPT  algorithm. 
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Figure  5.  The  basis  transformed  BTR-70  SAR  image  using  AWPT. 


4.2  Frequency-Dependent  Thresholding 


To  remove  clutters  in  the  AWPT  domain,  we  need  to  apply  a  thresh¬ 
old  level  to  the  transformed  coefficients  of  SAR  images.  The  common 
thresholding  method  is  the  soft-shrinkage  proposed  in  [10],  in  which  the 
above-threshold  signal  amplitude  is  reduced  by  a  level  of  the  thresh¬ 
old.  The  soft-thresholding  is  generally  applicable  to  the  cases  that 
signals  and  noise  or  clutter  are  spatially  overlapping.  For  SAR  im¬ 
ages,  the  target  image  and  clutters  are  basically  non-overlapping  even 
in  the  transform  domain.  Hard-thresholding  is  a  better  choice,  and  it 
is  defined  as  follows: 


[  5(m,n) 

if 

S(m,n ) 

1° 

if 

S(m,n ) 

>  T 
<  T 


(27) 


where:  T  is  the  threshold  level. 

If  the  spatial  clutters  are  white  and  Gaussian  distributed,  the  trans¬ 
formed  clutters,  as  we  have  demonstrated,  are  still  white  and  Gaus¬ 
sian  distributed.  Based  on  Newman-Pearson  Criterion  [15],  the  best 
threshold  level  is  a  constant  provided  that  it  meets  requirements  for 
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probabilities  of  false-alarm  and  target  detection.  But  for  SAR  images, 
a  pixel  has  some  weak  correlation  with  its  adjacent  pixels.  Therefore, 
the  clutter  on  a  SAR  image  is  not  really  white,  and  its  power  spectrum 
is  stronger  in  lower  frequencies.  We  design  the  threshold  levels  that 
are  dependent  on  the  central  frequency  for  the  quad-tree  decomposition 
outputs  at  each  of  the  final  branches.  For  the  output  of  one  of  out¬ 
ermost  branches  of  the  quad-tree  decomposition,  the  central  vertical 
and  horizontal  frequencies  are  assumed  to  be  fv  and  fjj,  respectively, 
and  the  central  frequency  at  that  output  branch  is  considered  to  be: 

fc  =  \JW+ftt  (28) 

From  quad-tree  decomposition  structure,  we  can  find  out  the  frequen¬ 
cies  fv  and  fa,  noting  that  due  to  the  down-sampling,  for  two- 
channel  decomposition  of  higher  frequency  band  data,  the  high  fre¬ 
quency  components  come  out  from  low  frequency  filter  output;  while 
the  low  frequency  components  from  high  frequency  filter  outputs. 
Therefore  the  threshold  level  for  the  quad-tree  decomposition  output 
branch  with  central  frequency  fc  is: 

™  =  (29) 

where  C  is  a  constant,  <jc  is  the  clutter  standard  deviation,  and  a 
and  j3  are  the  parameters  used  to  fit  the  mode  to  actual  clutters. 
Therefore,  the  threshold  level  is  the  largest  in  low  frequency,  and  the 
smallest  in  high  frequency.  Due  to  hard-thresholding  processing,  there 
is  no  DC  energy  loss  for  target  image  even  with  a  high  threshold  level 
in  near-DC  areas.  For  MSTAR  image  processing  we  choose:  a  — 
0.8  ~  1,  and  /3  =  0  in  (29)  to  approximate  the  clutter  spectrum 
characteristics  and  achieve  the  best  clutter  rejection  performance.  The 
AWPT  transformed  BTR-70  image  is  thresheld  with  parameters:  C  = 
0.5,  a  =  0.85,  and  f3  =  0  in  (29),  and  shown  in  Figure  6.  Using 
frequency-dependent  threshold,  we  still  can  keep  the  delicate  high- 
frequency  components,  which  represent  the  important  point-scattering 
reflections  in  the  SAR  images.  But  we  could  lose  or  blur  the  point¬ 
scattering  signals  after  the  processing  if  just  using  a  constant  threshold 
level  for  all  transformed  coefficients. 
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Figure  6.  The  transformed  BTR-70  SAR  image  after  frequency- 
dependent  thresholding  processing. 

4.3  Post-Processing 

Through  thresholding  processing  most  SAR  image  clutters  are  elim¬ 
inated  in  the  AWPT  domain.  The  original  SAR  images  can  be  re¬ 
stored  by  inverse-transforming  the  wavelet  coefficients  to  spatial  do¬ 
main  based  on  the  same  basis  in  the  quad-tree  decomposition.  Simi¬ 
larly  the  inverse-transform  can  be  efficiently  implemented  with  filter¬ 
ing  and  up-sampling  using  another  pair  of  filters  (P(n)}  and  (Q(n)}. 
Even  with  thresholding  processing  in  AWPT,  there  are  still  some  clut¬ 
ter  residues  existing.  After  the  inverse-transform,  those  clutter  resides 
will  spread  over  SAR  domain.  Therefore,  a  second  small  threshold  and 
a  clustering  processing  axe  applied  to  the  restored  SAR  image  to  fur¬ 
ther  improve  clutter  rejection  performance.  The  clustering  algorithm 
replaces  with  zeros  any  connected  non-zero  blocks  that  have  fewer  el¬ 
ements  than  a  pre-defined  block  size.  Usually  we  set  the  block  size  to 
larger  than  the  clutter  residue  sizes  and  smaller  than  target  block  size. 
For  MSTAR  images  a  clustering  block  size  of  32  is  usually  enough  to 
eliminate  all  clutter  residues.  Figure  7  is  the  restored  BTR-70  image 
through  post-processing  with  a  second  threshold  level  of  about  0.1<rc. 
Also  the  best  processing  result  using  direct  thresholding  and  clustering 
processing  is  shown  in  Figure  8.  Apparently  an  essential  part  of  the 
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Figure  7.  The  restored  BTR-70  image  with  thresholding  in  AWPT 
domain  and  post-processing  in  spatial  domain. 


Figure  8.  The  de-cluttered  BTR-70  image  with  direct  thresholding 
processing  in  spatial  domain. 
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target  is  lost  in  the  direct  thresholding  processing  because  of  low  signal- 
to-clutter  ratio  in  spatial  domain,  which  thus  poses  a  threat  to  correct 
target  identification  afterwards.  Through  the  AWPT  processing,  the 
whole  original  target  image  is  kept,  while  the  clutters  are  completely- 
removed. 

4.4  Performance  Metrics  for  MSTAR  Data  Processing 

For  a  SAR  image  recovery  processing  the  clutter  rejection  perfor¬ 
mance  of  an  algorithm  can  be  measured  based  on  the  ratio  of  mean 
square  error  to  clutter  standard  deviation.  But  in  the  cases  of  absence 
of  the  original  reference  signals  as  for  MSTAR  images,  the  improve¬ 
ment  of  average  target  image  signal-to-clutter  ratio  is  a  reasonable 
substitute.  We  define  average  target  Signal-to-Clutter  Ratio  (SCR) 
as: 

EE|s(«) 

SCR  =  '  3— -  (30) 

Nac 

where  S(i,j )  is  the  processed  target  pixel,  N  is  the  total  number  of 
target  image  pixels,  and  <rc  is  the  clutter  standard  deviation. 

The  clutter  rejection  performance  of  an  algorithm  can  be  mea¬ 
sured  based  on  SCR  improvements  through  its  processing.  In  the 
AWPT  clutter  rejection  processing,  we  gradually  increase  the  AWPT 
domain  threshold  level,  and  measure  the  SCR  improvements  due  to 
the  thresholding  processing  in  the  AWPT  domain  (prior  to  the  post¬ 
processing),  while  keeping  the  target  image  loss  inside  an  acceptable 
range.  Figure  9  shows  SCR  vs.  the  target  image  loss  for  three  different 
kinds  of  MTSAR  images  through  the  AWPT  processing.  Target  Image 
Loss  (TIL)  is  defined  as: 


TIL=!L=S  (31) 

where  5  and  S  are  the  average  target  signal  amplitudes  before  and 
after  the  processing,  respectively.  For  comparisons  the  SCR  vs.  TIL 
curves  for  Conventional  Wavelet  Transform  (CWT)  and  direct  thresh¬ 
olding  are  plotted  in  the  same  figures.  We  find  that  for  a  fixed  target 
image  loss,  the  AWPT  algorithm  almost  always  performs  better  than 
the  CWT  or  direct  thresholding  method. 
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Figure  9.  Signal-to-Clutter  Ratio  vs.  target  image  loss  using  AWPT, 
CWT  and  direct  thresholding  methods  for  (a)  BTR-70,  (b)  T-72,  and 
(c)  BMP-2  targets. 
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(C) 

Figure  9.  Continued. 

5.  CONCLUSIONS 

We  have  developed  a  new  clutter-rejection  algorithm  for  SAR  images 
based  on  Adaptive  Wavelet  Packet  Transform  (AWPT).  It  transforms 
the  basis  function  of  the  SAR  images  from  the  regular  pulse  basis  to 
a  wavelet  packet  basis  to  make  the  transformed  coefficients  maximally 
concentrated  on  the  transform  domain.  The  transformed  image  has 
a  higher  target  Signal-to-Clutter  Ratio  (SCR)  compared  with  that  in 
conventional  spatial  domain,  thus  the  clutters  can  be  removed  in  trans¬ 
form  domain  with  less  target  signal  losses.  The  method  is  based  on  the 
basic  assumption  that  a  clutter  pixel  in  a  SAR  image  is  not  or  weakly 
correlated  with  one  another;  while  target  signals  are  strongly  corre¬ 
lated  to  themselves.  Thus,  the  SCR  improvements  are  made  possible 
through  basis  transformation. 

The  energy  concentration  function  is  used  as  the  cost  function  to 
find  the  best  wavelet  packet  basis  to  make  the  transformed  coefficients 
with  the  maximum  SCR.  The  best  basis  search  and  the  basis  transfer- 
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mation  for  an  image  can  be  efficiently  implemented  using  2-D  quad-tree 
decomposition  algorithm.  To  search  the  best  wavelet  packet  basis  for 
a  SAR  image  with  clutters,  it  is  very  important  to  ensure  either  that 
the  clutters  and  the  target  image  are  approximately  non-overlapping 
after  transformation  in  transform  domain;  or  that  the  target  energy 
is  significantly  stronger  than  the  clutter  one.  Hence  the  best  wavelet 
packet  decomposition  trees  can  be  dominated  by  the  target  signals; 
otherwise  the  AWPT  algorithm  is  not  applicable. 

A  frequency-dependent  thresholding  method  is  introduced  because 
the  clutters  in  a  SAR  image  are  weakly  correlated,  i.e.,  they  are  “col¬ 
ored”  rather  than  “white,”  as  commonly  assumed.  The  threshold  level 
for  the  AWPT  transformed  outputs  are  chosen  based  on  the  central  fre¬ 
quency  of  the  data  in  this  output.  The  lower  the  central  frequency,  the 
higher  the  threshold  level.  Simulation  results  show  that  the  careless 
thresholding  could  degrade  the  performance  of  the  AWPT  algorithm 
to  that  of  conventional  wavelet  transform,  or  even  worse. 

Processing  results  on  MSTAR  images  show  that  this  algorithm  is 
very  effective  to  remove  the  clutters  for  a  SAR  image  with  only  lim¬ 
ited  clutter  information  such  as  standard  deviation  and  rough  spec¬ 
trum  characteristics  required.  With  robust  parameters  pre-selected, 
the  AWPT  algorithm  can  be  used  to  automatically  remove  the  clut¬ 
ters  in  a  large  number  of  SAR  images.  But  for  the  conventional  direct 
thresholding  method,  it  is  usually  very  difficult  to  find  a  threshold 
level  to  remove  all  background  clutters  and  keep  the  whole  target  im¬ 
age  unaltered.  This  new  method  is  also  very  useful  for  the  situations  in 
which  the  classical  direct  thresholding  method  might  not  work  at  all. 
Those  situations  include  that  some  parts  of  target  image  signals  are 
weaker  than  the  clutter,  or  that  the  target  signals  and  clutters  are  spa¬ 
tially  over-lapping  resulting  from  un-focused  processing,  or  multi-path 
effects. 
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Synthetic  Aperture  Radar  (SAR)  images  of  ground  targets  generally  consist 
of  target  features  clutters  from  background  scattering  [11-  In  automatic 
target  recognition  (ATR)  applications,  it  is  desirable  to  remove  the  clutter  from 
the  target  image  before  ATR  processing.  The  standard  way  to  suppress  chitteris 
to  apply  an  appropriate  threshold  level  to  the  whole  SAR  image-  However,  this 
uppffryrh  assumes  the  target  signaPto-dutter  ratio  (SCR)  is  large  enough. 
Otherwise  this  direct  threshold  approach  results  in  either  target  feature  loss  or 
remnant  clutter  residue.  In  this  work,  we  set  out  to  develop  a  dechxttering 
algorithm  to  automatically  extract  the  target  image  from  a  SAR  uny  by 
maximizing  the  SCR  using  the  adaptive  wavelet  packet  transform  (AWPT)  [2]. 
The  wavelet  N*<i<  is  a  generalization  of  the  conventional  wavelet  basis 

[3]  and  has  been  applied  fix  image  compression  [4]  and  moment  matrix 

sparsification  [5].  #. 

Our  approach  is  to  transform  the  SAR  image  to  a  new  domain  using  the 
wavelet  basis.  Sirce  a  typical  target  image  usually  consists  of  l»int 

scanners  more  diffused  region  features,  the  nnki-scaled  wavelet  basis  is 
well  suited  to  focus  the  target  image.  Clutter  image,  on  the  other  hand,  is 
statistically  uncorrelated  from  pixel  to  pixel,  and  the  transformed  clutter  nn&ge 
under  the  set  of  bases  remains  unfocused.  Therefore,  we  expect  that  the 
SCR  can  be  increased  by  transforming  the  original  image  by  an  appropriately 
chosen  set  of  wavelet  packet  basis.  The  cost  fiinction  of  our  AWPT  algorithm  is 
chosen  to  describe  bow  well  the  target  signal  is  focused  in  the  transform  domain. 
An  efficient  basis  search  algorithm  is  implemented  to  find  the  best  wavelet 
packet  basis.  Our  algorithm  is  tested  using  the  MSTAR  SAR  data  set  [6]  and 
results  show  th at  an  improved  SCR  can  be  achieved  using  die  AWPT  algorithm. 

2,  SAR  Image  Representation  with  Wavelet  Packet  Basis 
A  discrete  2-D  SAR  image  s(m,  n)  can  be  represented  as: 

r(m,n) « f(m,n)+c(m,n)  0 £n*m<N  0) 

where  t(m,  n)  and  c(m,  n)  denote  the  target  image  and  the  clutter  in  the  SAR 
image,  respectively.  We  define  a  set  of  orthogonal  and  complete  2-D  wavelet 
packet  basis  functions: 

0ij<J.0ip,q<2i.0ik.l<N2~1}  (2) 

where  j  denotes  the  scale  index.  J=logj(K).  p.  q  are  the  modulation  indices,  and 
k.  1  are  the  position  indices.  The  2-D  wavelet  packet  basis  function  can  be 
generated  from  the  product  between  two  1-D  wavelet  packet  bases: 
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y  is  a  l-D  wavelet  packet  basis  that  can  be  generated  from  the  scaling  function 
and  the  basic  wavelet  function  using  the  “2-scale  equation”  [3-5], 

The  transformation  of  a  SAR  image  sfm.  n)  using  (be  wavelet  packet  basis  is 
thus: 

$/,(*./)  =  2'n) 

m  m 

-SS»(m,«)£/^(i-2Vf-2M+2£c(«,n)C/^(i-2V/-2^) 

=f/,(*./)+C'|(i./) 

(4) 

where  T  ta&C  are  the  transform  coefficients  of  the  target  image  and  the  clutter 
m  the  image  with  the  wavelet  packet  basis. 

If  we  define  the  SCR  of  a  SAR  image  as  the  ratio  of  the  average  target 
amplitude  to  the  standard  deviation  of  the  clutter,  the  SCR  of  the  image  before 

and  after  the  transform  are  ?/<r(c)  and  f/a(C),  respectively.  The  clutter 
function  c (m.  n)  represents  the  reflectivity  of  different  pixels  in  the  background, 
and  &s  assumed  to  be  independent  and  identically  distributed  f7]- It  can  be  shown 
that  its  wavelet  packet  transform  coefficients  C  are  still  uncorrelated  [8]  and 
that  the  standard  deviation  of  the  clutter  does  not  change  after  the  basis 
tmsfonn.  On  the  contrary,  signals  from  the  target  area  are  strongly  correlated 
Md  itis  possible  to  increase  target  energy  concentration  with  a  basis 
teMMonnatjon.  Therefore,  with  a  given  SAR  image  we  need  to  search  and  find 
the  best  wavelet  packet  basis  to  maximize  the  amplitude  of  the  transformed 
target  image. 

3.  Adaptive  Search  Procedure  and  Implementation 

To  find  the  best  wavelet  packet  basis  and  implement  tbe  basis  transform,  we 
need  to  define  a  cost  function  to  describe  how  well  the  transformed  target  signal 
is  concentrated  with  a  wavelet  packet  basis.  The  best  basis  is  the  one  that 
generates  a  transformed  signal  having  tbe  least  cost  The  most  commonly  used 
cost  function  is  tbe  entropy  function.  Because  of  the  complexity  of  evaluating 
entropy  function,  we  use  the  energy  concentration  function  as  our  cost 
function  m  this  application.  For  a  transformed  SAR  image  it  is  defined  as: 

Co«=XZ^(i,/)|  (5) 

/•MU  1 

Because  there  are  many  possible  wavelet  packet  bases  in  (2)  for  the  transform,  it 
is  unpractical  to  tty  each  of  them  to  find  the  best  basis.  An  effective  quad-tree 
doomposmon  algorithm  was  generalized  from  that  proposed  in  [9]  to  find  the 
2-D  wavelet  packet  basis.  The  algorithm  decomposes  the  original  image 
using  2-channel  filtering  with  a  pair  of  quadrature  filters  through  all  scales  from 
foe  space  domain  to  foe  spectral  domain.  With  foe  cost  labeled  at  every  branch, 
foe  decomposition  tree  is  then  “pruned”  back  along  each  branch  from  the  last 
scale  toward  foe  earlier  scales.  The  pruning  process  is  accepted  whenever  it 
leads  to  a  lower  cost.  The  best  decomposition  tree  can  be  found  after  foe  search 
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process  with  a  total  computational  complexity  of  about  0(N2k>gN).  A  scale- 
dependent  threshold  is  then  applied  to  the  transformed  image.  Because  there  is  a 
weak  correlation  between  clatter  samples,  we  increase  the  threshold  level 
slightly  as  the  scale  increases.  The  threshold  level  is  chosen  as: 

Thj  =  Ka^jJJJ  (6) 

where  j  is  the  scale  index,  a  is  the  clutter  standard  deviation,  and  K  is  a  constant 
With  thresholding,  most  of  die  clutter  is  removed  in  the  wavelet  packet  basis 
domain,  and  the  image  is  inverse-transformed  back  to  the  SAR  domain  to 
restore  the  original  target  image  using  the  same  tree  from  the  decomposition. 
Although  the  SCR  is  much  improved  in  the  restored  image,  it  is  not  possible  to 
remove  the  clutter  completely  through  such  processing.  We  apply  a  very  small 
second  threshold  to  the  restored  image,  and  then  use  a  standard  clustering 
algorithm  to  get  rid  of  clutter  residues. 

4.  Test  Results 

To  test  the  effectiveness  of  the  AWPT  processing,  we  apply  it  to  the  MSTAR 
SAR  image  data  set.  Fig.  1  shows  an  MSTAR  image  in  which  the  target  is  a 
ground  vehicle  and  the  clutter  is  due  to  vegetation.  There  are  several  strong 
point  scatters  in  die  front  of  the  vehicle,  but  the  scattering  from  the  back  part  is 
relatively  weak.  Fig.2  shows  the  result  of  applying  the  direct  thresholding 
method  to  the  image,  fig.  3  shows  the  dechittered  image  by  applying  the 
AWPT  algorithm.  We  choose  Daubechies  filter  with  order  of  6  as  the  wavelet 
filter,  and  the  processing  parameters  used  are  K=1  and  0.05  c  for  the  second 
threshold.  By  visually  comparing  Figs.  2  and  3,  we  note  that  some  crucial 
features  of  the  target  are  kept  in  the  AWPT-processed  image.  In  both  processing 
methods,  there  is  some  target  information  loss.  Fig.  4  shows  the  signai-to-dutter 
ratio  versus  avenge  target  image  loss  for  the  two  processing  methods.  It  is 
observed  that  for  a  fixed  target  image  loss  the  AWPT  method  always  achieves  a 
higher  SCR  value  than  the  direct  thresholding  method.  Similar  results  are 
obtained  when  the  algorithm  is  applied  to  other  MSTAR  data. 
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Fig.  1 .  SAR  image  of  a  ground  Fig.  2.  Clutter  rejection  using 

vehicle  with  clutter  the  direct  thresholding  method 


Fig.3.  Clutter  rejection  using 
the  AWPT  approach. 


Fig.  4.  SCR  vs.  target  image  loss  for 
the  two  processing  methods. 
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Abstract 

Methodologies  for  simulating  the  radar  signatures  of  moving  ground  targets  using 
the  shooting  and  bouncing  ray  (SBR)  technique  ate  investigated.  The  motions  of  interest 
include  both  target  movement  as  a  whole  with  respect  to  background  and  sub-component 
motion  such  as  those  due  to  wheels  and  turrets.  The  brute-force  method  to  predict  the 
multi-pulse  radar  signature  is  to  simulate  the  target  range  profiles  at  N  different  dwell 
time  instances.  To  carry  out  the  brute-force  computation,  we  first  generate  N  CAD 
models  of  the  target  representing  the  states  of  the  target  and  compote  the  range  profiles 
for  all  N  models.  Not  only  is  the  first  step  a  tedious  process,  the  second  step  requires  a 
computation  time  that  is  N  times  longer  than  that  for  a  single  range  profile.  Furthermore, 
a  completely  new  run  must  be  carried  out  for  a  new  set  of  motion  parameters.  Therefore, 
two  alternate  approaches  are  investigated  to  more  rapidly  carry  out  dynamic  signature 
simulation. 

The  first  approach  is  based  on  an  extrapolation  algorithm.  It  is  an  extension  of  the 
signature  extrapolation  algorithm  for  target  look  that  we  have  developed  previously  for 
SBR  (Bhalla  and  Ling,  J.  Electromag.  Waves  Applications,  Feb.  1996).  In  this  approach 
we  launch  rays  only  once  at  the  target  Given  the  ray  history  and  the  target  motion,  we 
use  the  differential  Doppler  information  for  each  ray  to  predict  its  contribution  at  N  time 
instances.  The  extrapolation  algorithm  requires  only  a  single  ray  trace.  However,  a" 
completely  new  run  must  still  be  carried  out  for  a  new  set  of  motion  parameters. 

The  second  approach  is  to  use  the  extracted  3D  scattering  center  set  (Bhalla  and 
Ling,  IEEE  Trans.  Antennas  Propagat,  Nov.  1996)  from  the  target  to  predict  the  dynamic 
signature.  During  the  extraction  process,  each  scattering  center  is  tied  back  to  the. 
component  surfaces  on  the  target  that  gave  rise  to  that  scattering  center.  It  is  then 
possible  to  infer  the  motions  onto  the  individual  scattering  centers  associated  with  the 
component  Using  this  approach,  we  can  rapidly  predict  the  target  range  profiles  at  the  N 
time  instances.  In  addition,  it  is  easy  to  cany  out  this  procedure  for  arbitrary 
target/component  motions  at  little  additional  computation  cost  Dynamic  simulation 
results  for  ground  targets  with  various  motion  parameters  using  the  different  approaches 
will  be  presented. 
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ABSTRACT 

In  this  paper,  ISAR  images  generated  from  measured  data  are  compared  to  those  from  computer  simulation  in  order  to 
evaluate  the  effectiveness  of  ISAR-based  target  identification.  Three  sets  of  images  are  generated  including:  (i)  motion 
compensated  images  from  measured  data  using  a  joint  time-frequency  technique,  (ii)  reference  images  from  measured  data 
and  GPS-derived  aircraft  attitude  data,  and  (iii)  synthetic  images  predicted  by  Xpatch.  Visual  examination  and  correlation 
analysis  are  undertaken  to  compare  the  three  sets  of  images.  In  addition,  two  problem  areas  including  JEM  line  corruption  of 
the  measured  images  and  three-dimensional  rotation  of  the  target  are  identified. 

Keywords:  inverse  synthetic  aperture  radar  (ISAR)  imaging,  motion  compensation,  joint  time-frequency  technique,  Xpatch 
image  prediction 


1.  INTRODUCTION 

High-resolution  inverse  synthetic  aperture  radar  (ISAR)  imaging  has  been  regarded  as  a  possible  tool  for  target 
identification.1, 2  There  are  two  important  components  to  an  ISAR-based  target  ID  system.  The  first  component  is  the  image 
formation  algorithm  in  which  the  raw  data  collected  from  an  imaging  radar  is  processed  and  motion  compensated  to  form  a 
focused  image  of  the  unknown  target.  The  second  component  is  an  image  database  of  known  targets  populated  by  either 
actual  measurements  or  computer  simulation.  The  focused  image  from  measurement  is  then  matched  against  the  image 
database  in  order  to  determine  the  unknown  target  type.  The  success  of  an  ISAR-based  target  ID  system  is  therefore  critically 
dependent  on  the  quality  of  these  two  basic  building  blocks.  In  this  paper,  we  carry  out  a  comparison  of  ISAR  imageries 
generated  from  motion  compensated  measure  data  and  those  from  computer-simulated  synthetic  signatures.  Our  objectives 
are  to  provide  an  assessment  of  the  current  capabilities  and  identify  possible  hurdles  in  ISAR-based  target  ID. 

For  this  purpose,  three  sets  of  images  are  generated.  First,  the  motion  compensated  images  are  generated  from  radar 
measurement  data  of  an  airplane  in  flight.  Motion  compensation  is  carried  out  using  a  joint  time-frequency  technique  that  has 
been  reported  previously.4  We  shall  refer  to  these  images  as  the  JTF-mocomp  images.  Second,  a  set  of  reference  images  is 
generated  by  using  the  aircraft  motion  data  collected  during  the  flight  from  on-board  GPS  sensors.  Even  though  this  data  is 
available  only  from  cooperative  data  collection  and  not  in  the  real  target  ID  scenario,  it  serves  as  the  ground  truth  for 
evaluating  the  effectiveness  of  the  mocomp  algorithm.  We  shall  refer  to  these  images  as  the  GPS-reference  images.  Third, 
synthetic  images  are  simulated  from  a  CAD  model  of  the  aircraft  using  the  electromagnetic  signature  prediction  code 
Xpatch.5  We  shall  refer  to  these  images  as  the  Xpatch-synthetic  images. 

This  paper  is  organized  as  follows.  In  Section  2,  the  detailed  methods  are  given  on  generating  the  JTF-mocomp,  the  GPS- 
reference  and  the  Xpatch-synthetic  ISAR  images.  Results  are  presented  showing  example  images  and  the  correlation  between 
the  three  sets  of  images.  From  the  correlation  coefficients,  two  problem  areas  are  identified.  In  Section  3,  we  examine  the  jet 
engine  modulation  (JEM)  line  issue  in  the  frontal  look  region  of  the  target  and  propose  an  algorithm  to  remove  the  JEM  lines 
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from  the  measured  images.  In  Section  4,  we  examine  the  variable  imaging  plane  problem  during  certain  portions  of  the  flight. 
Conclusions  are  given  in  Section  5. 


2.  MEASURED  AND  SYNTHETIC  IMAGE  GENERATION 

With  the  radar  I/Q  data,  the  tracking  data,  the  GPS  data  and  the  CAD  model  of  the  aircraft,  measured  and  simulated  2-D 
ISAR  images  can  be  generated.  We  first  describe  the  procedures  used  to  generate  the  JTF-mocomp,  the  GPS-reference  and 
the  Xpatch-synthetic  ISAR  images.  From  these  results,  we  next  carry  out  a  comparison  to  evaluate  the  effectiveness  ofrthe 
motion  compensation  algorithm  and  the  electromagnetic  modeling. 

2.1.  Motion  Compensation  Using  Joint  Time-Frequency  Algorithm 

In  the  usual  case  of  ISAR  imaging  where  the  complex  target  motion  is  not  known,  motion  compensation  is  needed  to  form  a 
focused  ISAR  image.  For  this  purpose,  we  use  a  previously  developed  adaptive  JTF  algorithm.4  We  assume  that  after  the 
coarse  range  alignment,  all  the  scatterers  are  located  in  their  respective  range  cells.  The  radar  backscattered  signal  as  a 
function  of  dwell  time  tD  in  a  particular  range  cell  can  be  written  as 

E(tD)  =  Y  Ak  exp [~j—(R(tD)  +  xk  cos 0(tD)+yk  sin  0(tD))]  (l) 

"  c 

where  Nk  is  the  number  of  point  scatterers  in  that  range  cell,  and  Ak,  xk,  yk  are  respectively  the  scattering  amplitude,  down 
range  position  and  cross  range  position  of  the  k' *  point  scatterer.  R(tD)  is  the  residual  uncompensated  translation  displacement 
and  0(tD)  is  the  rotational  displacement.  The  JTF  technique  used  here  is  a  search  and  projection  procedure  to  represent  phase 
behavior  of  the  signal  E(tD).  To  find  the  motion  parameters,  basis  functions  in  the  form  of 

hit)  =  exp[-  j27t(f0t  +  /,f 2  +  J  f2t 3 )]  (2) 

are  chosen.  We  search  for  the  basis  function  over  the  parameter  space  (f0 ,  //,  /2)  that  best  represents  the  time-frequency 
behavior  of  the  signal  by  maximizing  the  projection  of  the  signal  onto  the  basis: 

max  \E(t)h*(t)dt  (3) 

foJhh 

After  the  time-varying  phase  for  the  strongest  point  scatterer  is  found,  we  multiply  the  original  signal  by  the  conjugate  of  this 
phase  factor  to  compensate  for  the  translation  motion.  This  algorithm  can  also  be  extended  to  multiple  range  ceils  to  correct 
for  higher  order  rotation  motion.  Figs.  1  and  2  show  the  JTF  processing  concept.  After  coarse  range  alignment,  a  particular 
range  bin  may  contain  multiple  scatterers  with  time-varying  Doppler  frequencies.  The  (dwell  time)-( Doppler  frequency) 
behaviors  of  these  scatterers  are  illustrated  in  Fig.  1(b).  Fig.  2(a)  shows  the  JTF  behavior  of  a  dominant  scatterer  in  a  range 
cell  from  actual  measurement  data.  We  see  from  Fig.  2(b)  that  after  the  JTF  compensation,  its  trajectory  is  straightened.  After 
applying  the  JTF  motion  compensation,  the  standard  FFT  processing  in  the  dwell  time  domain  brings  the  signal  into  the  cross 
range  image  domain. 

2.2.  Ground  Truth  Image  Generation  Using  GPS  Data 

The  motion  parameters  from  on-board  GPS  sensors  collected  during  the  cooperative  flight  of  the  airplane  are  next  used  to 
generate  the  ground  truth  ISAR  images.  The  resulting  reference  images  can  be  used  to  evaluate  the  quality  of  the  JTF- 
mocomp  images  from  the  last  section.  Furthermore,  the  motion  data  can  be  used  to  extract  azimuth  and  elevation  look  angle 
information  for  carrying  out  the  electromagnetic  simulation  of  the  aircraft. 

First,  coordinate  transformation  is  carried  out  since  the  GPS  data  is  in  the  fixed  Earth  system,  while  the  desired  azimuth  and 
elevation  angles  are  defined  with  respect  to  the  local  coordinate  of  the  aircraft  (see  Fig.  3).  The  latter  information  is  needed 
for  the  simulation  process  and  for  determining  the  absolute  scaling  of  the  measured  images  along  the  cross  range  dimension. 
Unlike  the  JTF  motion  compensation  where  the  coarse  range  alignment  is  carried  out  based  on  range  profile  correlation,  here 
range  data  from  the  GPS  measurement  are  used  directly  for  range  alignment.  In  addition,  the  aspect  angle  versus  dwell  time 


information  is  used  to  correct  for  higher-order  rotational  motions.  We  reformat  the  data  from  uniform  sampling  in  dwell  time 
to  uniform  sampling  in  aspect  angle.  The  FFT  is  then  used  to  generate  the  final  ISAR  images.  Since  the  look  angle 
information  is  available,  the  resulting  ISAR  images  can  be  scaled  in  the  cross  range  dimension  to  absolute  physical  size. 

2.3.  Synthetic  Signature  Prediction  Using  Xpatch 

To  test  the  effectiveness  of  electromagnetic  modeling  for  ISAR  imaging,  we  generate  the  simulated  ISAR  images  using 
Xpatch  and  the  radar  and  aircraft  motion  parameters  from  Section  2.2.  Xpatch  is  an  electromagnetic  computation  code  based 
on  the  shooting  and  bouncing  ray  method.6  It  can  be  used  to  compute  the  backscattering  of  complex  targets  of  large  electrical 
sizes.  In  the  Xpatch  simulation,  rays  are  shot  from  the  incident  look  angle  and  all  the  multiple  reflections  are  tracked  until  the 
rays  exit  the  target  (see  Fig.  4).  The  image  is  simulated  by  updating  the  ISAR  image  plane  one  ray  at  a  time  using  its  ray- 
spread  function.  The  image  update  is  further  accelerated  by  an  FFT-based  algorithm.  Note  that  this  fast  image  simulation 
algorithm  requires  only  a  single  ray  trace  per  image.  It  is  therefore  orders  of  magnitude  faster  than  the  conventional  method 
that  requires  multi-frequency,  multi-aspect  data.  The  typical  simulation  time  using  this  approach  is  approximately  30  minutes 
per  image  on  an  SGI  02  workstation. 

While  the  ISAR  images  generated  from  Xpatch  contain  absolute  scaling  along  the  cross  range  dimension,  the  ISAR  images 
from  the  JTF-mocomp  process  is  not  calibrated  along  that  dimension  since  the  rotational  speed  is  not  known.  Here  we  take 
advantage  of  the  look  angle  information  derived  from  the  GPS  data.  Prior  to  making  any  image  comparisons,  we  first  re-scale 
the  Xpatch  images  such  that  they  have  the  same  cross  range  resolution  and  scaling  as  the  images  from  the  measured  data. 

2.4.  Image  Comparisons 

In  this  section,  we  make  image  comparisons  among  the  three  sets  of  ISAR  images:  the  JTF-mocomp  images,  the  GPS- 
reference  images  and  the  Xpatch-synthetic  images.  Fig.  5  shows  one  set  of  such  comparisons.  The  look  angle  is  near  tail-on 
of  the  airplane.  The  JTF-mocomp  image  is  shown  in  Fig.  5(a).  It  is  fairly  well  focused  and  the  features  from  the  tail,  the  two 
wings  and  the  fuselage  are  clearly  exposed.  The  GPS-reference  image  is  shown  in  Fig.  5(b).  We  observe  that  the  agreement 
between  the  JTF-mocomp  and  the  GPS-reference  images  is  good.  Fig.  5(c)  shows  the  corresponding  Xpatch-synthetic 
image.  Again  the  outline  of  the  aircraft  is  readily  observed  in  this  image,  although  the  image  is  much  more  focused  than  the 
two  images  from  measured  data.  Our  experience  shows  that  the  JTF  motion  compensation  works  well  during  most  of  the 
fight  time.  After  a  side-by-side  visual  examination,  we  find  that  most  mocomp  images  are  in  fact  slightly  better  focused  than 
the  reference  images.  We  believe  this  is  due  to  the  limited  accuracy  of  the  GPS  sensor  data.  The  synthetic  images  show  good 
qualitative  agreement  with  the  measured  images.  However,  they  are  in  general  less  diffused  than  the  measured  images.  This 
is  not  surprising  since  the  synthetic  image  formation  assumes  no  motion  errors.  Furthermore,  the  CAD  model  used  does  not 
capture  all  of  the  fine  details  of  the  actual  target. 

While  the  comparison  from  visual  examination  shows  promising  agreement  between  the  three  sets  of  images,  an  image 
correlation  is  undertaken  for  a  more  quantitative  comparison.  Before  the  correlation,  the  images  are  power  transformed  to 
account  for  the  different  dynamic  ranges.  Fig.  6  shows  both  the  correlation  coefficient  between  the  JTF  and  GPS  images,  and 
that  between  the  JTF  and  Xpatch  images  versus  azimuth  look  angle.  From  the  two  curves,  we  see  that  the  JTF-mocomp 
images  agree  very  well  with  its  GPS-reference  counterparts,  indicating  that  blind  motion  compensation  is  a  very  feasible 
method  for  processing  real-world  radar  data.  The  correlation  coefficient  between  the  JTF  and  Xpatch  images  is  in  general 
lower  than  that  between  the  measured  images.  In  particular,  two  problem  regions  can  be  seen  from  this  plot.  First,  in  the 
region  near  nose-on  (180  degrees  in  azimuth),  the  correlation  coefficient  is  significantly  lower.  The  reason  is  due  to  JEM 
lines  in  the  measured  data.  This  problem  is  further  discussed  in  Section  3  and  an  algorithm  to  remove  JEM  lines  is  proposed. 
Second,  at  some  angles  around  the  broadside  region  (90  degrees  in  azimuth),  the  correlation  coefficient  is  also  low.  The 
associated  JTF  images  are  found  to  be  of  low  quality.  After  further  investigation,  it  is  found  that  the  image  blurring  is  due  to 
variations  in  the  imaging  plane,  not  the  motion  compensation  processing.  This  problem  is  discussed  in  detail  in  Section  4. 


3.  JEM  REMOVAL 

Jet  engine  modulation  is  a  phenomenon  caused  by  the  high-speed  rotational  movement  of  the  aircraft  engine.9,  10  For  an 
imaging  radar,  the  typical  PRF  is  much  slower  than  the  engine  rotation  frequency.  Thus  the  resulting  ISAR  image  in  the 
frontal  region  of  an  aircraft  contains  an  aliased  component  along  the  cross  range  dimension.  Such  effect  is  difficult  to  predict 
accurately  using  simulation.  Furthermore,  JEM  lines  are  noise-like  and  can  corrupt  the  geometrical  features  of  the  target  in 


the  ISAR  image.  For  target  ID  using  2D  ISAR  imaging,  it  would  be  useful  to  devise  an  algorithm  to  remove  JEM  lines,  and 
therefore  enhance  the  image  and  the  subsequent  ID  process. 

Let  us  assume  that  the  aircraft  consists  of  a  slow  moving  body  with  a  constant  rotational  velocity  and  a  fast  moving 
engine  component  with  a  different  rotational  velocity  Then  the  received  radar  return  as  a  function  of  dwell  time  can  be 
written  as: 

A  4 

E(tD)  =y.Ak  exp [-j - ( R(tD)  +  xk  cos (QbtD)  +  yk  sin(£VD)] 

tri  c 

(4) 

N  A'trf  V  7 

+  X  A*  exP[-^'  —(r(*d  )  +  xk  cos( aptD )  +  yk  sin(Q ptD )] 

k=Nb+ 1  C 

where  N  is  the  total  number  of  point  scatterers  within  one  range  cell,  of  which  Nb  are  the  body  scatterers.  Usually  Qp  is  much 
greater  than  Qb.  While  the  first  term  can  be  meaningfully  mapped  into  the  image  plane  of  the  target  via  the  Fourier  transform, 
the  second  term  results  in  serious  Doppler  smearing  across  the  cross  range  domain  and  may  overshadow  the  target  features. 

We  can  also  utilize  the  joint  time-frequency  technique  described  earlier  under  motion  compensation  to  separate  the  fast 
moving  part  from  the  relatively  slow  moving  body.11  For  the  component  due  to  target  body  scattering,  the  Doppler  frequency 
is 


fo  =—  0A[ycos(Qfc?D)  +  xsin(QArD)]*^QA(y  +  xQferD)  (5) 

c  c 

while  the  Doppler  frequency  due  to  the  fast  rotating  part  is 

fo  =  —  Q„[;ycos(Q  f0)  +  xsin(Q  tD)]  (6) 

c 

We  can  see  that  (5)  is  a  linear  function  of  dwell  time  while  (6)  is  a  sinusoidal  function.  In  the  time-frequency  plane,  the  two 
signals  can  thus  be  distinguished.  If  we  further  parameterize  the  signal  by  basis  functions  that  have  linear  Doppler  frequency 
behavior  as  a  function  of  dwell  time,  the  two  signals  can  be  separated  automatically  by  their  displacement  and  slope 
parameters.  We  utilize  the  adaptive  joint  time-frequency  processing  technique  to  carry  out  the  parameterization.  The  basis 
used  is  similar  to  that  given  in  (2)  with  the  linear  and  quadratic  phase  terms.  The  project  and  search  procedure  given  in  (3)  is 
carried  out  iteratively.  At  each  iteration,  the  basis  parameters  (f0>  //)  and  Bp,  which  is  the  maximum  projection  value  of  the 
signal  onto  the  basis,  are  found.  The  best  basis  at  stage  p  is  then  removed  from  the  signal: 

Ep+l(t)  =  Ep(t)-Bphp(t)  (7) 

The  searching  process  is  iterated  until  the  energy  of  the  residue  signal  is  smaller  than  a  preset  threshold.  The  signal 
component  due  to  the  target  scattering  can  thus  be  reconstructed  by  using  all  the  bases  with  small  displacement  f0  and  small 
slope  parameter  fr.  Fig.  7  shows  the  correlation  between  the  synthetic  images  and  the  measure  images  after  JEM  removal.  We 
observe  that  the  correlation  coefficients  in  the  JEM  region  are  increased  after  we  remove  the  JEM  interference  from  the  body. 


4.  IMAGING  PLANE  VARIATION 

From  the  correlation  curves,  some  images  generated  from  the  measured  data  outside  of  the  JEM  angular  region  are  also  found 
to  be  poorly  focused.  For  those  frames,  we  have  found  that  the  associated  GPS-reference  images  also  fail.  If  we  assume  the 
data  quality  from  the  radar  did  not  change  abruptly,  the  cause  must  not  be  in  the  motion  compensation  algorithm,  but  rather 
from  extraneous  motions  during  the  imaging  interval,  which  cannot  be  handled  by  any  motion  compensation  method.  Here 
we  examine  the  possible  causes.  The  movement  of  the  aircraft  in  space  relative  to  the  ground  radar  consists  of  radial  motion 
and  rotational  motion.  Because  the  range  alignment  process  can  remove  the  radial  motion,  what  remains  a  problem  is  the  3D 
rotational  motion  of  the  aircraft.  Suppose  during  the  imaging  interval  the  rotational  motion  is  described  by  the  look  angles  (0, 
0)  on  the  target.  The  requirement  of  a  constant  rotation  axis  described  by  (0O,  <t>o)  leads  to  the  following  equation  that 
constraints  the  values  of  (0,  <()): 


(8) 


sin  0O  sin  6  cos  (<j>  -  0O )  +  cos  cos  6  =  0 

If  0o=O,  eq.  (8)  requires  0=n/2.  In  this  case,  the  aircraft  rotates  about  its  vertical  axis  and  the  imaging  plane  is  the  top  view  of 
the  airplane.  If  0o=n/2  and  <|>o=7c/2,  then  <j)=0.  In  this  case,  the  aircraft  rotates  about  its  axis  along  the  wing  and  the  imaging 
plane  is  the  side  view  of  the  airplane. 

The  images  from  these  two  cases  will  be  quite  different.  If  an  imaging  interval  includes  both  types  of  motion,  the  resulting 
image  is  expected  to  be  unfocused.  Although  this  is  an  extreme  example,  we  do  observe  several  cases  in  the  real  data  where 
the  imaging  plane  variation  in  the  128  pulse  records  leads  to  bad  images.  Fig.  8  shows  a  particular  example  of  imaging  plane 
variation.  From  the  azimuth -elevation  plot  in  Fig.  8(a)  derived  from  the  GPS  data,  we  recognize  two  different  imaging 
planes.  During  the  first  half  of  the  128  records,  the  aircraft  undertakes  nearly  a  vertical  rotation,  while  during  the  second  part 
it  undertakes  a  horizontal  rotation.  So  we  get  a  very  smeared  image  as  shown  in  Fig.  8(b).  We  can  get  a  better  image  by  using 
only  the  first  part  or  the  second  part  of  the  data.  The  resulting  image  from  the  first  part  of  the  data  is  a  side  view  of  the 
aircraft,  while  the  image  from  the  second  part  of  the  data  is  a  top  view  of  the  aircraft.  Unfortunately,  we  usually  do  not  have 
access  to  the  aircraft  attitude  data  on  non-cooperative  targets.  In  such  cases,  the  questions  of  how  to  detect  imaging  plane 
variation  and  how  to  form  the  best  possible  images  become  important  research  issues.  This  topic  is  currently  under 
investigation  by  us. 


5.  CONCLUSIONS 

To  evaluate  the  effectiveness  of  ISAR-based  target  ID,  ISAR  images  generated  from  measured  data  are  compared  to  those 
from  computer  simulation.  Three  sets  of  images  are  generated  including:  (i)  motion  compensated  images  from  measured  data 
using  a  joint  time-frequency  technique,  (ii)  reference  images  from  measured  data  and  GPS-derived  aircraft  attitude  data,  and 
(iii)  synthetic  images  predicted  by  Xpatch.  Visual  examination  and  correlation  analysis  are  undertaken  to  compare  the  three 
sets  of  images.  Through  the  comparisons,  the  following  conclusions  can  be  made.  First,  JTF  motion  compensation  method 
performed  well  with  real  world  ISAR  data.  Second,  Xpatch  is  a  feasible  tool  to  generate  a  synthetic  image  database  for 
ISAR-based  target  ID.  Two  problem  areas  are  also  identified.  JEM  line  corruption  of  the  measured  images  is  quite  severe  in 
the  frontal  sector  of  air  targets.  A  possible  algorithm  based  on  joint  time-frequency  technique  is  proposed  to  remove  the  JEM 
lines.  In  addition,  the  problem  of  imaging  plane  variation  is  identified  to  be  the  cause  of  poor  ISAR  images  during  some  time 
of  the  flight.  Further  investigation  is  needed  to  devise  methods  to  overcome  this  problem. 
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Fig.  3.  Global  and  local  coordination  systems  of  the 
problem  (figures  not  drawn  to  scale). 


Fig.  4.  Shooting  and  bouncing  ray  technique  for  radar 
signature  simulation. 
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Fig.  8.  Effect  of  imaging  plane  variation  on  1SAR  image. 


Publication  [13] 


Application  of  Adaptive  Joint  Time-Frequency 
Processing  to  ISAR  Image  Formation 


Hao  Ling  and  Junfei  Li 
Dept,  of  Electrical  and  Computer  Engineering 
The  Univ.  of  Texas  at  Austin 
Austin,  TX  78712-1084 


I.  Introduction 

High-resolution  inverse  synthetic  aperture  radar  (ISAR)  imaging  is  a  promising  tool  for  non-cooperative 
target  identification  (NCTI).  The  main  challenge  in  ISAR-based  NCTI  is  to  form  a  well-focused  image  of 
an  articulating  target  with  unknown  motion.  In  this  paper,  we  first  review  the  application  of  joint  time- 
frequency  methods  for  ISAR  image  formation.  By  using  an  adaptive  joint  time-frequency  (AJTF)  algorithm 
to  estimate  the  phase  of  the  prominent  scatterers,  we  show  that  the  target  motion  can  be  estimated  and  a 
focused  image  of  the  target  can  be  constructed.  Results  of  applying  the  algorithm  to  measured  ISAR  data 
are  presented  and  discussed.  Secondly,  we  report  on  our  recent  work  to  extend  the  AJTF  algorithm  to 
address  the  more  challenging  situation  when  the  motion  of  the  target  is  not  limited  to  a  two-dimensional 
plane.  In  particular,  we  discuss  our  research  to:  (i)  detect  the  presence  of  three-dimensional  motion  using 
the  AJTF  algorithm,  and  (ii)  develop  algorithms  to  focus  the  image  in  the  presence  of  complex  three- 
dimensional  motion. 


2.  ISAR  Motion  Compensation  Using  Joint  Time-Frequency  Algorithm 

We  first  review  the  application  of  joint  time-frequency  methods  for  ISAR  image  formation.  To  form  a 
focused  image  from  raw  radar  data,  it  is  customary  to  first  carry  out  a  coarse  alignment  of  the  data  in  the 
range  dimension,  followed  by  fine  motion  compensation  in  the  cross  range  dimension.  Joint  time- 
frequency  techniques  have  been  shown  to  be  a  useful  tool  to  carry  out  the  fine  motion  compensation  [1,2]. 
We  assume  that  after  the  coarse  range  alignment,  all  the  scatterers  are  located  in  their  respective  range 
cells.  The  radar  backscattered  signal  as  a  function  of  dwell  time  t  in  a  particular  range  cell  can  be  written  as 

N  Arrf 

E(t)=  £  Akexp[-j—^-(R(t)  +  xkcosd(t)  +  yksinO(t))]  (1) 

k=l  c 

where  N  is  the  number  of  point  scatterers  in  that  range  cell,  and  Ah  xkl  yk  are  respectively  the  scattering 
amplitude,  down  range  position  and  cross  range  position  of  the  k*  point  scatterer.  R(t)  is  the  residual 
uncompensated  translation  displacement  and  Oft )  is  the  rotational  displacement.  Due  to  translation  and 
rotational  motion,  the  Doppler  frequency  versus  dwell  time  behavior  of  the  point  scatterers  within  this 
range  cell  is  not  constant  in  the  joint  time-frequency  plane  (see  Fig.  1).  An  effective  JTF  technique  to 
extract  the  motion  parameters  is  based  on  a  search  and  projection  procedure  to  represent  the  phase  behavior 
of  the  signal  Eft).  This  procedure  is  based  on  the  adaptive  spectrogram  proposed  in  [3],  and  is  similar  in 
concept  to  a  one-term  matching  pursuit  algorithm  [4].  We  shall  term  it  the  adaptive  JTF  (AJTF)  algorithm. 
To  find  the  motion  parameters,  basis  functions  in  the  form  of 

h(t)  =  exp[-j(ajt  +  a2t2jra3t3)]  (2) 

are  chosen.  We  search  for  the  basis  function  over  the  parameter  space  (aIy  a2>  a3)  that  best  represents  the 
time-frequency  behavior  of  the  signal  by  maximizing  the  projection  of  the  signal  onto  the  basis: 


max  |j  E( t  )h* ( t  )dt 

dj  tQj 


2 


(3) 


After  the  time-varying  phase  for  the  strongest  point  scatterer  is  found,  we  multiply  the  original  signal  by 
the  conjugate  of  this  phase  factor  to  compensate  for  the  translation  motion.  This  algorithm  can  also  be 
extended  to  multiple  range  cells  to  correct  for  higher  order  rotation  motion.  After  applying  the  JTF  motion 
compensation,  the  standard  FFT  processing  in  the  dwell  time  domain  brings  the  signal  into  the  cross  range 
image  domain.  Results  of  applying  the  algorithm  to  simulated  and  measured  ISAR  data  will  be  presented 
and  discussed. 


Range  profiles  versus  pulse  no. 
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Fig.  1.  Fine  motion  compensation  is  carried  out  by  extracting 
the  Doppler  frequency  versus  dwell  time  behavior  of 
the  strong  point  scatterer  in  the  signal. 


3.  Three-Dimensional  Motion  Estimation  Using  Joint  Time-Frequency  Algorithm 

One  basic  assumption  of  standard  motion  compensation  algorithms  is  that  the  target  only  undergoes 
motion  in  a  two-dimensional  plane  during  the  dwell  duration  needed  to  form  an  image.  From  several 
independent  examinations  of  measured  ISAR  data  sets  recently,  it  was  reported  that  the  presence  of  three- 
dimensional  motion  is  quite  detrimental  to  focusing  the  image  [5-7].  We  shall  report  on  our  recent  work  to 
extend  the  AJTF  algorithm  to  address  the  more  challenging  situation  when  the  motion  of  the  target  is  not 
limited  to  a  two-dimensional  plane.  In  particular,  we  discuss  our  research  to:  (i)  detect  the  presence  of 
three-dimensional  motion  using  the  AJTF  algorithm,  and  (ii)  develop  algorithms  to  focus  the  image  in  the 
presence  of  complex  three-dimensional  motion. 
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Fig.  2.  (a)  Simulated  2D  target  motion,  (b)  Phase  behavior  of  the  prominent  point 
scatterer  in  range  cell  1  extracted  using  AJTF.  (c)  Phase  behavior  of  the  prominent 
point  scatterer  in  range  cell  2  extracted  using  AJTF.  (d)  Ratios  of  the  extracted 
phase  parameters  from  the  two  range  cells.  Note  that  they  are  nearly  constant. 
(e)-(h)  Similar  to  (a)-(d),  except  that  3D  motion  is  assumed.  The  resulting  ratios  in 
(h)  are  no  longer  constant. 


Allowing  for  arbitrary  three-dimensional  motion  in  space,  we  consider  the  following  model  as  a 
generalization  of  the  model  for  two-dimensional  motion  in  ( 1 ): 

N  47lf 

E(t)=  exp[-j-l±(xk  +  ykO  +  zk</>)  1  (4) 

k=l  c 

where  0is  the  azimuth  angle  of  the  target  with  respect  to  the  radar,  and  <p  is  the  elevation  angle.  In  (4),  it  is 
assumed  that  the  translation  motion  has  been  removed  and  that  the  standard  small-angle,  small  bandwidth 


approximations  apply.  This  model  reduces  to  the  standard  two-dimensional  motion  model  when  0  and  <f> 
are  linearly  related. 

In  general,  a  focused  image  cannot  be  obtained  from  the  standard  two-dimensional  motion 
compensation  algorithm  when  three-dimensional  target  motion  is  present  due  to  model  mismatch. 
Therefore,  it  would  be  useful  to  detect  the  presence  of  three-dimensional  motion  directly  from  the  radar 
data.  Our  approach  is  to  utilize  the  AJTF  algorithm  to  extract  the  phase  behavior  of  the  radar  data  at 
multiple  range  cells.  We  first  parameterize  the  phase  of  the  prominent  point  scatter  in  one  range  cell  using 
AJTF.  Next  we  repeat  the  same  procedure  at  another  range  cell.  It  can  be  shown  that  when  the  target 
undergoes  only  two-dimensional  motion  during  the  dwell  duration,  the  ratio  between  the  parameters  (ah  a2, 
a3)  extracted  from  one  range  cell  and  those  corresponding  parameters  in  another  range  cell  should  be 
constant.  Therefore,  by  examining  the  ratio  of  the  parameters,  we  can  distinguish  two-dimensional  motion 
from  three-dimensional  motion.  Fig.  2  illustrates  the  idea  using  simulated  point  scatterer  data.  Figs.  2(a)- 
(d)  show  the  two-dimensional  motion  scenario  and  Figs.  2(e)-(h)  show  the  three-dimensional  scenario.  It 
can  be  seen  from  the  results  in  Fig.  2(d)  that  the  determined  ratios: 

Ci  =at{ range  cell  l)/at{ range  cell  2)  (5 ) 

are  nearly  constant  for  all  the  terms  in  case  of  two-dimensional  motion,  as  expected.  For  three-dimensional 
motion,  the  ratios  are  not  the  same,  as  seen  in  Fig.  2(h).  Results  of  applying  this  concept  to  radar  data  in 
order  to  isolate  the  time  durations  when  the  target  undergoes  three-dimensional  motion  will  be  presented. 
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